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ABSTRAGT 

Turbulence is ubiquitous in molecular clouds (MGs), but its origin is still unclear because MGs 
are usually assumed to live longer than the turbulence dissipation time. Interstellar medium (ISM) 
turbulence is likely driven by SN explosions, but it has never been demonstrated that SN explosions 
can establish and maintain a turbulent cascade inside MGs consistent with the observations. In 
this work, we carry out a simulation of SN-driven turbulence in a volume of (250 pc)^, specifically 
designed to test if SN driving alone can be responsible for the observed turbulence inside MGs. We 
find that SN driving establishes a velocity scaling consistent with the usual scaling laws of supersonic 
turbulence, suggesting that previous idealized simulations of MG turbulence, driven with a random, 
large-scale volume force, were correctly adopted as appropriate models for MG turbulence, despite the 
artificial driving. We also find that the same scaling laws extend to the interior of MGs, and that the 
velocity-size relation of the MGs selected from our simulation is consistent with that of MGs from the 
Outer-Galaxy Survey, the largest MG sample available. The mass-size relation and the mass and size 
probability distributions also compare successfully with those of the Outer Galaxy Survey. Einally, 
we show that MG turbulence is super-Alfvenic with respect to both the mean and rms magnetic- 
field strength. We conclude that MG structure and dynamics are the natural result of SN-driven 
turbulence. 

Subject headings: ISM: kinematics and dynamics - MHD - stars: formation - turbulence 


I. INTRODUCTION 

Understanding molecular cloud (MG) turbulence is key 
to understanding the star formation process, because su¬ 
personic turbulence is ubiquitous in MGs and drives their 
fragmentation into stars. Supersonic turbulence has been 
studied extensively in the context of star formation, and 
its statistical properties are at the core of recent mod- 


els of the star formation rate dKrumholz & McKee||2005 

Padoan & Nordlund||20IIl 

Hennebelle & Chabrier||2Ull 

Eederrath & Klessen||20I2 

) and the stellar initial mass 

function ( 

Padoan et al.||l9y7l Padoan & Nordlund||2002 

Hennebel] 

e Chabrier||2UU8 Hopkins 2U12p. Most nu- 


merical studies of supersonic turbulence have used a ran¬ 
dom large-scale force to drive the turbulence, as custom¬ 
ary in the turbulence literature. It remains to be shown 
that such an idealized external force is a good approx¬ 
imation of the actual large-scale processes driving the 
interstellar-medium (ISM) turbulence. While the simu¬ 
lations use a volume force that penetrates the interior of 
the fluid, as long-range forces do (e.g. gravity and mag¬ 
netic fields), the real ISM driving forces may be surface 
forces, such as large-scale shocks from spiral arms or su¬ 


pernova (SN) bubbles. The effect of different types of 
surface forces, or of a combination of volume and surface 
forces on the turbulence and, thus, on star formation, 
has not been studied systematically. 

The turbulence in MGs is also key to understanding 
their origin. The generation and maintenance of MG tur¬ 
bulence must be an integral part of the cloud formation 
process, because most of the energy of an MG is in the 
form of turbulent kinetic energy. Eor example, the great 
majority of small and intermediate-mass MGs are known 
to have rather larg e virial parameters (see §10 and Heyer 


et al. 


2001, 2009[), thus their turbulent energy is large 


n to form and disperse them in a few dynamical 
The same may be true also for the most massive 


enougn 
times. 

MGs, even if their virial parameter tends to be of order 
unity. The spatial and velocity structures of MGs fol¬ 
low power laws that span all scales from the smallest to 
the most massive clouds, suggesting a universal origin of 
clouds and cloud turbulence. 

In this work, we adopt the viewpoint that MG turbu¬ 
lence is just one component of the general ISM turbu¬ 
lence. The question of the origin of MG turbulence is 
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thus turned into the more general question of how the 
ISM turbulent energy is shared among its different gas 
phases. If we demonstrate that the total large-scale tur¬ 
bulent energy of the ISM is dynamically consistent with 
the turbulent energy in its dense phase, no extra energy 
source specific to MCs is needed. 

It is generally accepted that SN explosions dominate 
the energy budget of star-forming galaxies at MC scales, 
although l arge-scale gravitational instabilities in galactic 
disks (e.g. Elmegreen et al.||2QQ3 Bournaud et al.||2QlQ|) 


and gas compres sion m spiral density waves (e.g. 5e- 
.|2015|) may also contribute to the turbulence, 
mt-bchmi 


menov et al. 


'I'he Kennicut-bchmidt star formation law of disk galax¬ 
ies, despite giving a large gas-consumption timescale of 
order 1 Gyr, corresponds to an energy input from SN 
explosions that exceeds the turbulence dissipation rate 
of the ISM of those galaxies. The analysis of HI maps 
of nearby face-on galaxies also leads to the conclusion 
that SN feedback is responsible fo r the observed HI line 
width, except in the disk outskirts ([Tamburro e t al.|2QQ9; 
Stilp et al.|2Q13 lanjamasimanana et al. |2015| . Detailed 


modeling of the disk vertical balance has also shown that 
SN feedback can maintain the ISM turbulence that de¬ 
termines the dis k scale height, resulting in self-regulated 
star formatio n dOstriker et a h |2Q1 Q ©striker & Shetty 
^iguere et al.|| a picture that may 


iguere et al 




2011 Faucher 

apply also to gala xies with high star f ormation rate at 
redshift z = 1 — 3 (jLehnert et al.||2013|) . 

Prior to these studies, galactic-fountain simulations 
had already demonstrated that SN explosions can drive 


the observed ISM turbulen ce (de Avillez & Breitschw- 
erdt|2005 Joung et al.|2009). However, these simulations 
barely resolved the formation of the most massive clouds. 


and could not resolve their internal structure and dynam¬ 
ics in order to compare with observed MC properties. 
Furthermore, the evolution of individual SN explosions 
was not resolved with high-enough spatial resolution to 
study their interaction with individual MCs. 

Recent numerical works have studied the momentum 
injection by individual SN explosions into the ambient 
medium, assumin g realistic ISM density and tempera¬ 
ture fluctuations (|Walch fc Naab 2Q15| Martizzi et al. 


2015 Iffrig & Hennebellej 2015 "Kim & ©striker 2U15p 
'I'hese studies are useful for the derivation of feedback 


models to be implemented as sub-grid physics in galaxy 
formation simulations, but do not address the problem 
of the generation and maintenance of MC turbulence by 
SN explosions. 

In this work, we focus on a smaller-scale region than 
in galactic-fountain simulations, and use high-enough nu¬ 
merical resolution to model the MCs, their internal struc¬ 
ture and dynamics, and the evolution of individual SN 
bubbles with sub-pc resolution (see §2). We also select 
clouds from the simulation, as connected regions above 
a threshold density (§3), to study their turbulence and 
to compare with observed MCs. First we study statis¬ 
tical properties over the whole computational volume, 
such as total energies (§4), power spectra (§5) and ve¬ 
locity structure (§6.1). Then we present properties of 
individual clouds selected from the simulation, such as 
the velocity structure (§6.2), the virial parameter (§7), 
the cloud lifetime (§8) and the magnetic field strength 
(§9). Finally, we compare the properties of the clouds 
selected from the simulation with those of MCs from the 


Five College Radio Astronomy ©bservatory (FCRA©) 
©uter Galaxy Survey (§10). 

2. NUMERICAL SETUP 

Th e simulation is carried out with the Ramses A MR 
code (|Teyssier|2002|). We refer to|Padoan et al.|(|2014|) for 
a brief description of our version of Ramses. What is n ew 
here, relative to the simulation injPadoan et al. (2014), is 
the use of the full energy equation (instead of assuming 
an isothermal equation of state), and the inclusion of SN 
explosions and tracer particles, besides the much larger 
physical size of the computational volume. 

We simulate a cubic region of size Lbox = 250 pc (large 
enough to contain a few turbulence correlation lengths, 
while small enough to allow sub-pc resolution), with a 
mean density of 5 cm“^ (corresponding to a column den¬ 
sity of 30 MqPc”^ and a total mass of 1.9 x 10^ Mq) and 
a Type H SN rate of 6.25 Myr“^ (or a galactic rate of 
100 Myr“^kpc“^, if all SN explosions occurred within the 
vertical extent of our box). We do not consider Type la 
SNe, because of their lower rate and higher scale height, 
and because we distribute SN explosions randomly, so 
our SN rate could also be interpreted as the sum of the 
Type H and Type la rates. 

These rate and column density values are consistent 
with the Kennicutt-Schmidt relation, and our compu¬ 
tational volume may be viewed as a dense section of a 
spiral arm. For example, the total column dens ity in the 
Perseus arm o f the Milky Way is 23 MqPc”^ (Heyer & 


Terebey 1998). However, we do not include a galactic 


gravitational field, and adopt periodic boundary condi¬ 
tions in all directions, so vertical stratification and out¬ 
flows of hot gas are neglected. We have chosen this ide¬ 
alized setup because one of the motivations of this work 
is to relate the statistical properties of SN-driven turbu¬ 
lence to previous studies of randomly-driven, supersonic, 
isothermal turbulence that were carried out on periodic 
boxes without stratification. Furthermore, we relate the 
velocity scaling to theoretical predictions that also ne¬ 
glect complications such as gravity and stratification^ 
Besides the pdV work, and the thermal energy in¬ 
troduced to model SN explosions, our en ergy equation 
adopt s uniform photoelectric heating as in |Wolfire et al. 
(1995), with efficiency e = 0.05 and the FUV radiation 


held of Habing (1968) with coefficient Go = 0.6, cho¬ 
sen to obtain temperature distributions co nsistent with 
those from the comprehensive simulations by |Walch et al.| 
(2015). Because the code conserves total energy, kinetic 
and magnetic dissipations are included self-consistently 
as energy sources (this dissipation is purely numerical, 
as we do not include viscosity or resistivity explicitly). 
We use a tabulated optically thin cooling fu nction con¬ 
structed frorn the extensive compilation by [Gnedin 

(|2012|), based on 7 5 million runs of the Cloudy 


Hollon 


code (|Ferland et al.||1998| to sample a large range of 
conditions, and from whicn the results have been made 


publicly available as a Fortran code with accompanying 
database. All relevant atomic transitions are included in 
the Cloudy runs. Although available in Cloudy, molecu- 


^ Star-formation simulations with gravity have shown that, un¬ 
der reasonable conditions (e.g. MCs not c ollapsing as a whole). 

gravi ty does not affect the velocity scaling (jFederrath & Klessenj 

- 
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Fig. 1. — Phase diagram of gas pressure versus density based 
on 11 snapshots covering uniformly the time period when gravity 
is included in the simulations, t = 45 to 56 Myr. The gray scale 
represents the square root of the volume fraction at a given pressure 
and density. 

lar cooling is not included because the runs are restricted 
to a single computational zone, with a negligible column 
density, to enforce the optically thin case. Above a tem¬ 
perature of 100 K, atomi c cooling is dominant up to den¬ 
sities of 10^ cm“^ (e.g. Neufeld et al. 1995). At lower 
temperature and high densities, molecular cooling should 
be included. However, molecular cooling and cosmic-ray 
heating are neglected, and their thermal balance in very 
dense gas is emulated by clamping the resulting drop in 
tempe rature at 10 K. As pointed out by [Gnedin fc Hollon 


(2012), including the balance of molecular cooling and 


cosmic-ray heating in such a treatment does not make 
much physical sense, since the balance of these processes 
at high densities crucially depends on radiative trans¬ 
fer effects; in particular the absorption of UV radiation 
by small-scale high-density cloud structures. In the ab¬ 
sence of radiative transfer, and given the optically thin 
assumption, we approximate the UV shielding in MCs 
by tapering off the photoelectric heating exponentially 
above a number density of 200 cm“^ (assuming a char¬ 
acteristic size of 1 pc for MC structures at our critical 
density, c orresponding to a cr itical visual extinction of 
0.3 mag MFranco fc Cox[l986| ). 

Figure shows the phase diagram of gas pressure ver¬ 
sus density sampled over the last 11 Myr of the simu¬ 
lation. The horizontal feature around densities of a few 
10“^^ gcm“^ is a consequence of the approximation of 
self-shielding. In a real MC, or a model where the ab¬ 
sorption of UV-radiation by the filamentary structure of 
dense gas and dust is taken into account more realis¬ 
tically, the transition to opaqueness would take place at 
different densities at different locations, and similar local 
phase diagram features would be washed out . Comparing 
our cu rrent phase diagram with the ones in jWalch et al.j 
(2015|), we actually find the largest discrepancy not to be 
near that horizontal feature, but rather at higher densi¬ 
ties where, for some reason, their balance of heating and 
cooling results in a temperature of about 30 K, instead 
of our assumed value of 10 K. 


The simulation is started with zero velocity, a uniform 
density nH,o = 5 cm“^, a uniform magnetic field Bq = 
4.6 fiG and a uniform temperature Tq = 10^ K. The 
first few SN explosions rapidly bring the mean thermal, 
magnetic and kinetic energy to approximately steady- 
state values, with the magnetic field amplified to an rms 
value of 7.2 jaG. The value of the mean magnetic field 
is chosen to achieve near equipartition with the kinetic 
energy at large scales, as shown in §4. It also yields an 
average value of \B\ of 6.0 /iC, consistent with the value 
of 6.0 ± 1.8 juG derived from th e ‘Millennium Arecibo 21 - 
cm Absorption-Line Survey’ by Heiles & Troland (2005). 

The simulation is run for 56 Myr (with a total of 359 
SN explosions), initially without tracer particles and self¬ 
gravity. At t = 33 Myr we include 150 million passively- 
advected tracer particles following the mass distribution 
in the computational volume (each particle represents 
a fluid element of 0.013 Mq). Tracer particles are ad- 
vected with the same symplectic Kick-Drift-Kick scheme 
used for dark matter particles in Ramses, but the kicks 
are ignored - they are passive - and velocities in the drift 
are instead sampled using CIC interpolation of fluid ve¬ 
locities. 

At t = 45 Myr we include self-gravity. Several clumps 
of order 1,000 Mq start to collapse. Larger resolution 
and sink particles are required at that point, which will 
be the topic of a followup work. For the purpose of this 
work, we stop the simulation after approximately 11 Myr 
of evolution with self-gravity. This is also motivated by 
the need to trace the exact location of SN explosions 
when clouds have been forming massive stars for a time 
of order 10 Myr, as explained in the next section. 

To enforce a sub-pc resolution of the evolution of SN 
bubbles and of their interaction with the ISM, we adopt 
AMR criteria based on density, density gradients and 
pressure gradients. Although our root grid contains only 
128^ cells, our AMR criteria result in rather large volume 
filling factors of high-resolution cells: 75% of our com¬ 
putational domain is covered at a resolution equivalent 
to 256^, 22% at a resolution of 512^, and 2% at 1024^. 
Because the volume filling factor of the clouds selected 
in this work is approximately 0.5%, our clouds are all 
resolved at the maximum resolution. 

The left panel of Figure shows the logarithm of the 
projected density of the whme computational volume, at 
the final time of the simulation. The structure is highly 
filamentary and appears to be self-similar. The right 
panel of Figure shows a sub-region magnified by a fac¬ 
tor of four; the same type of filamentary structure is seen 
at this smaller scale. This structure is very similar to that 
previously found in supersonic simulations of isothermal 
turbulence, and consistent with the appearance of nearby 
GMCs mapped with the Herschel satellite. 

Apart from SN explosions, the simulation neglects any 
other energy source, such as winds and radiation feed¬ 
back from massive stars, or external forces, such as spi¬ 
ral arm shocks and fluctuations of the large-scale gravi¬ 
tational potential, that may affect MC turbulence in real 
galaxies. This choice allows us to test if SN turbulence 
alone can explain the origin and maintenance of CMC 
turbulence, while also providing a significant reduction 
of the computational cost. We also neglect the chem¬ 
ical evolution of the ISM. Although we do not expect 
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Fig. 2.— Left Panel: Logarithm of projected density along the x-axis of the simulation volume. The mean magnetic field direction is 
horizontal on the image. The column density value is larger in darker regions, with black set to a maximum value of 5 x 10^^ cm“^, and 
white to a minimum value of 6 x 10^® cm“^. The actual values of column density span a much larger range, and these limits have been 
chosen to optimize the contrast of the density structure. The image includes the whole computational volume, that is a size of 250 pc. 
Right Panel: Same as in the left panel, but for a region four times smaller (62.5 pc) shown by the white squared in the left panel. The 
color-table limit values have been increased by a factor of 4 and 5, 2 x 10^^ cm“^ (black) and 3 x 10^^ cm“^ (white), to match the larger 
mean density in this region. 


that the dynamics would be strongly affected by a more 
precise computation of cooling and heating based on dy¬ 
namically evolved chemical abundances, the selection of 
MCs and their comparison with observations would cer¬ 
tainly benefit from a dynamical computation of the H2 
and C O abundances (jGlover & ClarkJ[2012| |Walch et al. 


times, neglecting the possibility of spatial and tempo- 
ral clustering. In the galactic fo unta in simulations by|d^ 


2014). The neglect of chemistry in this work is solely 


motivated by considerations of computational cost. 

2.1. SN Explosions 

Individual SN explosions are implemented with an in¬ 
stantaneous addition of 10^^ erg of thermal energy and 15 
Mq of gas, distributed according to an exponential pro¬ 
file on a spherical region of radius rsN = ^dx = 0.73 pc, 
wh ere dx is our smallest cel l size, dx = I/box/1024 = 0.24 
pc. ‘ ' 

tioh 


Kim & Ostriker (2015) have derived a useful condi- 
tor numerical convergence of the evolution of SN 
remnants, which states that both the grid size and the 
initial SN radius must be smaller than one third of the 
shell-formation radius, dx^rsN < rsf/3. Because i n our 
case rsN > dx, an d given the expression for rgf in |Kim 
& Ostriker (2015), the condition for our simulation is 
where no is the ambient density. 


FSN 


< lOpcnQ 


Given our value of tsn = 0.73 pc, the condition becomes 
no ^ 300 cm“^. The volume filling factor of gas at den¬ 
sity above that value varies in the approximate range be¬ 
tween 0.0003 and 0.0009, with the largest value reached 
only at the end of the simulation, while the total num¬ 
ber of SN explosions in the simulation is 359. Because 
the locations of SN explosions are randomly distributed, 
there is only a small probability that one or more SN ex¬ 
plosions violate the condition for numerical convergence. 

SN explosions are generated at random positions and 


Avillez & Breitschwerdt (2005) and Joung et al. (2009) 
it was assumed that 60% of the SN explosions occurred 
in clu sters, and only the remaining 40% at random loca¬ 
tions. Joung et al. (2009) assumed that the clustered SN 
population is distributea in clusters with up to 40 stars 
more massive than 8 Mq, A^sn = 40, with a probability 
of cluster size following a power law, ^ A'g^. All massive 
stars from a given cluster were assumed to explode at the 
same location, though at different times distributed in an 
interval of 40 Myr, the approximate lifetimes of the least 
massive stars to explode. 

However, assuming that clusters have a typical velocity 
dispersion of the cold gas of order 10 km/s, and consid¬ 
ering the average value of A^sn = 10, the typical cluster 
would contribute on average 1 SN explosion every 4 Myr, 
with a separation of approximately 40 pc between con¬ 
secutive SN expl osions, due to the cluster motion not 
accounted for by Joung et al. (2009). The 10 explosions 
would cover a distance of 400 pc over 40 Myr . Thus, even 
assuming the reasonable cluster statistics of [Joung et al.| 
(20091), the explosions would not appear to be clustered. 
A realistic prediction of spatial and temporal clustering 
of SN explosions requires that the formation and kine¬ 
matics of individual stars is numerically resolved through 
the adoption of very high spatial resolution and sink par¬ 
ticles. In the absence of that, a random SN distribution 
is a reasonable assumption. 

Another potentially important issue is the correlation 
between the location of SN explosions and the position 
of their parent clouds. Recent simulations have shown 
that the large-scale structure of the ISM is sensitive to 
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the amount of correlation between SN positions and den¬ 


sity p e aks |G atto et al.|2Q15| |Walch et al.||2014|) . |Walch 
et al. ( 2Q14[ ) have concluded that the random dikribu- 


tion case or the case of clustering in space and time, but 
not correlated with density pea ks, are favored by the ob¬ 
serv ations. On th e contrary, in Dobbs & Pringle (2013) 
and Dobbs (|2015 ) all SN explosions are assumed to oc¬ 
cur inside MCs. As an MC (or any density peak above a 
certain threshold) is created, SN feedback is turned on in 
its interior. This unrealistic SN feedback does not allow 
to address the question of the origin of MC turbulence, 
a s this is driven directly by SNe by design. 


Iffrig & Hennebelle (2015) have addressed this issue by 
focusing on the ettect of a single SN on a single MC of 
10^ Mq, showing that the effect of a SN explosion on a 
MC is much stronger if the explosion occurs inside the 
cloud than outside of it. Their results may be affected by 
their specific choice of the ambient density, 1.2, 20 and 
700 cm“^, for the three SN positions they tested. The 
probability of high density is certainly increased inside 
a MC than outside of it. However, the probability of a 
SN explosion at high density must be quite small, be¬ 
cause the volume filling fraction of dense filaments and 
clumps is low even within the volume enclosing a MC. 
Furthermore, HII regions and winds would probably pre¬ 
vent SN explosions in dense gas in general. Nevertheless, 
it seems reasonable that SN remnants expanding from 
within a MC would be more effective at bringing mate¬ 
rial above the escape velocity than SN remnants pushing 
on a MC from the outsidn^ _ 

With a Salpeter IMF (|Salpeter|1955), the mean stellar 
mass between 8 and 100 Mq is approximately 19 Mq, 
corresponding to a stellar lifetime of approximately 9 
Myr. Assuming that it takes at least 1 Myr to initi¬ 
ate star formation in a young cloud, and to accrete the 
stellar mass, the SN feedback would thus be important 
in a MC after a characteristic time of 10 Myr (assum¬ 
ing that massive stars can remain in the general region 
of their parent cloud for a parent-cloud crossing time). 
We will show in §8 that our most massive clouds of or¬ 
der 10^ Mq selected towards the end of the simulation 
were formed approximately 20 Myr earlier (tform ~ 20 
Myr in Figure [21]). For these clouds, thus, the SN feed¬ 
back from locally-formed massive stars should start to 
play an important role, depending on the delay between 
cloud formation and the formation of the first massive 
stars. 

The same is true for all clouds with lifetime of order of 
10 Myr or larger. Because we find in §8 that the cloud 
lifetime is on average four times longer than the cloud 
dynamical time, and using the expression (25) for the 
dynamical time derived in §10.4, the condition is that 
the dynamical time is at least 5 Myr or longer, or, equiv¬ 
alently, that the cloud mass is larger than 500 Mq. As 
the formation of massive stars most likely requires clouds 
more massive than 500 Mq as well, we conclude that SN 
feedback from locally-formed massive stars should gen¬ 
erally play a role for most clouds forming massive stars. 
As mentioned in §2, this is in fact the main reason why 
the simulation was stopped approximately 11 Myr after 
introducing self-gravity, as neglecting the locally-formed 
massive stars beyond that time would be unrealistic. The 
effect of the correlation of SN explosions with their parent 
clouds will be considered in a future work where the for¬ 


mation and kinematics of individual massive stars will be 
numerically resolved in order to model self-consistently 
the precise time and location of SN explosions. 

3. MC SELECTION 

The main question addressed by this work, whether 
SN explosions can drive and maintain the observed MC 
turbulence, can be partly answered independently of the 
definition of MCs, by computing the velocity structure 
functions for the dense and cold gas throughout the com¬ 
putational volume. As long as most of such gas is in 
clouds, its global structure functions should be equiva¬ 
lent to the average of the structure functions of MCs. 
Nevertheless, it is important to characterize the turbu¬ 
lence within individual clouds and to also compute other 
cloud properties that can be compared with observations. 

Because chemistry is not included in this work, we can 
only define MCs as cold over-densities in the SN-driven 
ISM turbulence. A detailed comparison with the obser¬ 
vations is beyond the scope of this work. It would require 
synthetic observations, hence molecular abundances from 
chemical-network calculations. In the absence of syn¬ 
thetic observations, we prefer to avoid the selection of 
MCs in position-position-velocity (PPV) space, and in¬ 
stead define MCs in 3D space (PPP) in the simplest pos¬ 
sible way, as connected regions above a single threshold 
gas density, nH,min- In order to test if our results de¬ 
pend on spatial resolution or threshold density, we ex¬ 
tract clouds using three different mesh sizes, 2dx, Adx 
and Sdx (we create uniform grids of 512^, 256^ and 128^ 
cells, respectively), and four different threshold densi¬ 
ties, nH,min = 100, 200, 400 and 800 cm“^, generating 
12 cloud catalogs. Examples of clouds from one such cat¬ 
alog (nH,min = 100 cm“^ and 512^ resolution) are shown 
in Figure The images show the projected density of 
18 clouds, the 2nd to the 7th most massive ones in each 
of three snapshots. 

The analysis of MC properties derived with the full 
three-dimensional information, such as the results on MC 
velocity scaling, the discussion on the virial parameter 
and cloud structure, the evaluation of cloud lifetimes and 
the study of the cloud magnetic field, is based on clouds 
extracted with nH,min = 100 cm“^ and 128^ resolution. 
For this analysis, the resolution only affects the definition 
of cloud boundaries, as all results are derived from the 
position, velocity, density and magnetic field of the tracer 
particles, thus taking advantage of the highest resolution. 

When we compare with observational data, deriving 
projected quantities such as surface density, equivalent 
radius and line-of-sight velocity dispersion, we verify the 
results on all 12 catalogs. For the mass and size distri¬ 
butions, we also report quantitatively on the dependence 
of the slope of their power-law tails on cloud extrac¬ 
tion density and resolution. However, all the plots we 
show in the comparison with the observations are based 
on the highest-resolution catalog (512^ cells, or 0.49 pc) 
and on the threshold density that best matches the ob¬ 
served mass-size relation, nH,min = 200 cm“^. Further¬ 
more, velocity dispersions are based on tracer particles 
and so take advantage of the highest available resolution, 
dx = 0.24 pc0 This resolution matches well the high- 

^ The number density of tracer particles is very large in dense 
gas, where the mesh is refined to the highest resolution. 
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Fig. 3. — Square root of the projected density of 18 clouds extracted with nH,min = 100 cm“^ and 512^ resolution from three different 
snapshots near the beginning, the middle, and the end of the time interval with self-gravity. The clouds are numbered in order of decreasing 
mass. To save space in the panels, the most massive cloud from each snapshot (cloud 0) is not shown. The grey color table covers a range 
of column densities from 0 Mqpc“^ (black) to 200 Mqpc“^ (white). The column density of gas bel ow t he thre shold density of 100 cm“^ 
is not included. Cloud 6 of the top panel and cloud 4 of the bottom one are the same as in Figures and [T^ respectively. The complex 
filamentary structure is very similar to that observed in MCs. 
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Fig. 4.— Total kinetic, magnetic and thermal energies in the 
simulation versus time, during the period analyzed in this work, 
between t = 33.06 Myr and t = 56.43 Myr. The energies are mea¬ 
sured in 188 snapshots (8 snapshots per Myr), integrating over the 
whole simulation volume, independent of density. The horizontal 
dotted line shows the initial magnetic energy (also the magnetic 
energy corresponding to the mean magnetic field in the simula¬ 
tion, which is a conserved quantity). There is near equipartition 
between the lowest values of kinetic energy, the mean value of the 
magnetic energy, and the highest values of the thermal energy. 


est one in the observational survey we consider (Heyer 
et al.||2QQl ). After selecting clouds with circular velocity 


Vc < 20 km s“^ and mass Md > 100 Mq, we are left 
with 3,228 observed clouds with measured distances cor¬ 
responding to a range of spatial resolutions of 0.24 — 3.0 


pc. 

We select from the simulations only clouds more mas¬ 
sive than 100 Mq to guarantee that, even at the lowest 
value of nH,min = 100 cm“^, the smallest clouds contain 
more than 1,000 computational cells (our largest cloud 
of nearly 3 x 10^ Mq contains more than 3x10^ cells). 
Because we use 150 million tracer particles, our small¬ 
est and lowest-threshold-density clouds contain a mini¬ 
mum of approximately 7,000 particles, while our most 
massive cloud of 3 x 10^ Mq contains more than 20 mil¬ 
lion particles. With this minimum mass, each simulation 
snapshot yields over 200 clouds. In the comparison with 
the observations (see §10) we use 7 snapshots from ap¬ 
proximately the last 6 Myr of the simulation, to include 
the effect of self-gravity, resulting in sample sizes ranging 
from 595 clouds in the smallest catalog (128^ resolution 
and nH,min = 800 cm“^), to 1,615 in the largest one 
(512^ resolution and nH,min = 100 cm“^). The plots we 
show in §10, based on 512^ resolution and nH,min = 200 
cm“^, use measurements from 1,547 clouds. The catalog 
sizes could be considered three times larger, as projected 
quantities are computed in the three orthogonal direc¬ 
tions. However, the clouds (and cloud masses) would be 
the same, so the three samples would not be completely 
independent. Thus, we compare with the observations 
using only one of the directions, after verifying that the 
results are independent of direction. 

Despite the mass range of over three orders of mag¬ 
nitude of our clouds, we refer to all of them as MCs, 
instead of following the common nomenclature that clas¬ 
sifies clouds as cores, clumps, molecular clouds, giant 
molecular clouds, and giant molecular cloud complexes. 


Fig. 5. — Same as Figure^ but with the energies computed only 
for gas with density above lOO cm“^, the same threshold density 
adopted for the cloud selection. There is a clear energy separation 
in the dense gas, with the lowest kinetic energy values being ap¬ 
proximately an order of magnitude larger than the mean magnetic 
energy, and two orders of magnitude larger than the lowest values 
of the thermal energy. Thus, the turbulence in the dense gas is 
both supersonic and super-Alfvenic. 

in order of increasing mass. As we view all clouds as 
cold density enhancements of the ISM turbulence, and 
because of the scale-free nature of the turbulence, that 
nomenclature is not useful for this work. 

4. TOTAL ISM ENERGIES 

We analyze our simulation in a time interval start¬ 
ing after the turbulence has been fully developed until 
the end of the simulation, between t = 33.06 Myr and 
t = 56.43 Myr. This is also the time interval for which 
we have tracer particles in the simulation. The evolu¬ 
tion of the total kinetic, magnetic and thermal energies 
in that time interval is shown in Figure The three 
energies are all around 10^^ erg. Both the thermal (Eth) 
and kinetic (E'k) energies show strong oscillations, while 
the magnetic energy (E'm) is almost constant with time. 
The kinetic energy oscillates mostly above 10^^ erg, and, 
interestingly, at its valleys, the magnetic energy is al¬ 
most in equipartition with E\^. The thermal energy is 
the smallest, mostly below both the magnetic and ki¬ 
netic energies. 

Given the strong temporal oscillations in kinetic en¬ 
ergy, it may seem surprising that the magnetic energy 
does not experience significant fluctuations at all. The 
reason is that the kinetic energy peaks correspond to the 
energy injected by SN explosions mainly in the form of 
expansions and shocks. While the magnetic field can be 
strongly compressed and amplified by the shocks, the vol¬ 
ume filling factor of the postshock gas is tiny, and thus 
the total magnetic energy of the computational volume is 
barely changed, as illustrated by the following estimate. 

Consider a SN remnant of radius i?sN and a preshock 
magnetic field Bi. Compression by the SN shock would 
amplify the magnetic field strength by a factor equal to 
the density-jump factor, T, of the shock. In the adia¬ 
batic phase, T = 4 and in the radiative phase, T AIa 
with At A = Vsh/{Bi/ ^y4:7^pl) the Alfvenic Mach num¬ 
ber of the shock. As the remnant expands, solenoidal 
turbulent motions also develop (see §5.1), which can am¬ 
plify the magnetic field as well. The timescale for the 
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magnetic energy amplification by solenoidal motions is 
roughly the turnover time of the largest edd ies in the 
postshock region (e.g. Federrath et al. 2011a). We as¬ 
sume that the turbulent velocity in the postshock region 
is the post shock velocity, Ugh/F, and the large eddy size 
is the thickness of the post shock region of the remnant, 
which is i?sN/(3r) (due to the compression by the 
shock). The large eddy turnover time is then estimated 
to be i?SN/(3ush)- Considering the age of the remnant 
is ^ |i?SN/'^sh and fi^sN/'^sh in the adiabatic and ra¬ 
diative phases, respectively, solenoidal motions may am¬ 
plify the magnetic energy by one e-fold or so, meaning 
an amplification factor of A 3 for the magnetic en¬ 
ergy. Including both the effects of shocks and solenoidal 
motions, the magnetic energy density in the post shock 
region would be -^AV‘^B‘1. 

The above description of the magnetic field amplifi¬ 
cation only applies to the compressed layer behind the 
shock. Due to the compression by a factor of F, the width 
of the compressed layer is given by i?sN/(3r), so the total 
magnetic energy within the compressed layer is estimated 
to be On the other hand, the magnetic en¬ 

ergy in the hot cavity interior to the compressed layer is 
small due to the expansion and can be neglected. Consid¬ 
ering that the magnetic strength is still Bi in regions 
not reached by the SN remnant, the total magnetic en¬ 
ergy in the simulation box is -^BKAttATNs^R^^/S^ 

Uox - 47rA^sN-RiN/3), where A^sn is the number of SNe 
exploded around the same time and Wox is the volume of 
the simulation box. Because the highest kinetic-energy 
peaks are of the order 10^^ erg, A^sn may be as large as 
^10. If we define a filling factor of each SN remnant 
as / = 47ri?gjy^/(3kbox), the total magnetic energy can 
be written as ^[A^sN(^r — 1)/ + IJ^iWox, meaning 


that is amplified by a factor of A^sN(^r — 1)/ + 1. 
In the Sedov phase, F is constant, and the amplification 
factor for the magnetic energy increases with the filling 
factor, /. Applying the physical conditions in our simu¬ 
lation, we find that, when the Sedov phase ends due to 
radiative cooling, / increases to 3 x 10“^. Thus, with 
F = 4 and A 3^ (AF — I)/ is only 0.003. Therefore, 
due to the tiny filling factor of the SN remnant, the total 
magnetic energy is amplified only by a negligible amount 
during its early evolution. 

When the remnant evolution enters the radiative 
phase, the jumping factor, F ^ AIa, which can be 
significantly larger than 4. Using the pressure-driven 
snow-plough solution to account for the deceleration of 
the shock velocity (and hence the decrease of Ma with 
time) and the increase of i?sN (and the filling factor, 
/), we find that the maximum of the amplification fac¬ 
tor, (A^F — I)/, by a single SN in the radiative phase is 
0.05. Therefore, even if at a given time there are 10 
SNRs that happen to simultaneously amplify the mag¬ 
netic energy by the maximum possible amount, the total 
amplification of the magnetic energy is only 50%. Note 
that, in the above estimate, we have ignored the dynam¬ 
ical effect of the magnetic field and the dissipation of the 
magnetic energy, and thus the realistic amplification due 
to the SN explosions may be considerably smaller than 
50%. This explains the near constancy of the total mag¬ 
netic energy despite the strong oscillations in the total 
kinetic energy of the flow. 



Fig. 6.— Velocity dispersion versus time. The rms velocity is 
computed over the whole computational volume (dashed line) or 
only for gas with density larger than 100 cm“^ (solid line). The 
dotted line shows the mass-weighted rms velocity, averaged over 
the whole volume. 


Although the SN energy is introduced purely as ther¬ 
mal energy in the simulation, the thermal energy is effi¬ 
ciently converted into kinetic energy, and is also partly 
radiated, so the resulting turbulence is mildly supersonic, 
with a time-averaged kinetic to thermal energy ratio of 
(F’k/^th) =3.75. A crucial question for this work is what 
fraction of that total kinetic energy is given to the dense 
gas and, thus, is available as turbulent kinetic energy 
for the MC turbulence. The answer is given in Figure 
showing the time evolution of the total energies in the gas 
with density larger than 100 cm“^. The time-averaged 
ratio of dense-gas kinetic energy and total kinetic energy 
is (F^k,d/^k) = O.II. The ratio of the kinetic energies per 
unit mass is ((F^k,d/Afd)/(F^k/Afbox)) = 0.55 (the dense- 
gas mass fraction is (Md/Mbox) = 0.18), and oscillates 
in time between 0.19 and 0.97, with the largest values 
achieved during the peaks of total kinetic energy. This 
large ratio means that, per unit mass, the kinetic energy 
is only a factor of ^ 2 less than equally distributed be¬ 
tween the dense gas and the rest of the gas. This high ef¬ 
ficiency of kinetic energy transfer to the dense gas results 
in a realistic velocity dispersion of dense gas, (crv,d) = 8.5 
km/s, as shown in Figure]^ 

Figure also shows a clear separation of energies in the 
dense gas, with (F;k,d/^m,d) = 27.4 and (F;m,d/^th,d) = 
9.7. Thus, the turbulence in the dense gas is both super¬ 
sonic and super-Alfvenic, as further confirmed below for 
individual clouds sele cted from our simulation (see ^10 ) 
and first suggested by Padoan & Nordlund (1997 1999). 


5. POWER SPECTRA AND DRIVING SCALE 

The expansion of SN remnants deposit energy on a 
broad range of scales, affecting the velocity scaling of 
the turbulence. We observe this in the form of wave-like 
features in the velocity power spectrum or deviations in 
the velocity structure functions (see §6) at small and in¬ 
termediate scales in nearly one fourth of the snapshots of 
our simulation, as expected for our rate of approximately 
6 SN explosions per Myr, and a characteristic expansion 
time to 20-40 pc of ^ 4 x 10^ yr. 

Figure shows the power spectra from a single snap¬ 
shot that captures the early expansion phase of a single 
SN remnant (here it has reached a diameter of approx- 
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Fig. 7.— Power spectra of velocity, E(k), solenoidal velocity 
component, Es{k), compressive velocity component, Ec{k), square 
root of kinetic energy, E}<i{k), and magnetic field, Em(k)^ obtained 
from the average of the power spectra of the x, y and 2 : components 
of these fields in the root grid (128^ computational cells) of a single 
snapshot at time t = 52.75 Myr. This time captures the early 
expansion of a SN remnant to a diameter of approximately 30 pc, 
as indicated by the peak of the velocity power spectrum. The wavy 
appearance of Ec{k) is explained in the text. The fact that Eq does 
not show similar fiuctuations, and E^ < Ec at the energy peak, 
indicates that this remnant has so far expanded into a relatively 
uniform hot medium, where the baroclinic effect is negligible and 
solenoidal modes are not efficiently generated yet. 


imately 30 pc) that has not had any major interaction 
with dense gas yet, so most of the freshly-injected energy 
is still in the compressive modes (see §5.1)|^ The Figure 
shows the power of the velocity, E{k)^ of the solenoidal 
(divergence-free) velocity component from a Helmholtz 
decomposition, Es{k), of the compressive (curl-free) ve¬ 
locity component from the same Helmholtz decomposi¬ 


tion, Ec{k), of the square root of kinetic energy (p^/^r^i), 
EK{k), and of the magnetic field, Ejji{k). They are com¬ 
puted from the root grid (128^ computational cells) of 
a single snapshot at time t = 52.75 Myr. The velocity 
power spectra, especially the compressive one, exhibit 
significant fiuctuations at large and intermediate k. The 
wavy behavior is a signature of strong SN shocks and 
can be understood with a one-dimensional illustration. 
Consider two SN shocks moving in opposite directions 
away from the origin (the explosion center). If the ve¬ 
locity profile in between the two shocks is assumed to 
be roughly linear with the distance from the origin, it 
is straightforward to show that the power spectrum is 
oc k~‘^{cos{kR) — sm{kR)/kR)‘^^ where R is the radius 
of the SN shock. Clearly, this spectrum oscillates at 
k > R~^, which explains the wavy behavior of the veloc¬ 
ity spectra shown in Figure 
In between SN explosions, or in regions not directly 
affected by a rapidly expanding SN remnant, the flow 
should have time to relax from the transient state with 
direct SN impact, and we expect the velocity spectra 
to be smoother. To test this, we consider a time range 
t = 40.5 — 48.8 Myr, around the middle of our integration 
time. This choice avoids the larger SN rate towards the 


Fig. 8. — Power spectra as in Figure ^ but averaged over 66 
snapshots within the time interval t = 40.5^48.8 Myr. During this 
time interval the random SN rate is a bit lower than during the first 
and last third of the simulation, so there is a reduced probability 
that a snapshot captures the early expansion of SN remnants with 
the corresponding perturbations to the power spectra. As a result, 
the time average is representative of the statistically relaxed power 
spectra. The power-law slopes are obtained from a least-square fit 
in the range of wave numbers 4.2 < fcLbox/27r < 12.1. 

beginning and the end of the simulation (see Figure]^, 
and thus minimizes the number of snapshots with direct 
impact from very recent SN explosions. Figure shows 
the average three-dimensional power spectra computed 
from the root grid of 66 snapshots in the chosen time 
range. For each snapshot, the 3D spectra are obtained 
by averaging the three components (e.g., the average of 
the power spectra of Ux^ Uy and Uz). The average spec¬ 
tra we obtain are qualitatively similar to those of the 
most relaxed snapshots, and they exhibit a clear (though 
short) inertial range. The slopes have been computed in 
the wavenumber interval 4.2 < /cI/box/27r < 12.1, corre¬ 
sponding to the scale interval 59.3pc > 20.6pc, where 

all spectra are very well described by power laws. The 
measured inertial-range slopes are given in Figure To 
our knowledge, this is the first time that inertial-range 
slopes are identified in SN-driven turbulence. Because 
of the importance of this result, a full discussion of the 
power spectra will be presented elsewhere. Here, we fo¬ 
cus only on two points of direct interest for this work, the 
ratio of compressive to solenoidal modes and the driving 
scale. 

5.1. Ratio of Compressive to Solenoidal Modes 

SN driving should not be viewed as a purely com¬ 
pressive form of driving. Our simulation indicates that, 
around the effective driving scale (see §5.2), compres¬ 
sive modes are not dominant over solenoidal ones. We 
find that the compressive-to-solenoidal ratio is typically 
smaller than unit y, as shown by the time-averaged power 
spectra in Figure [Sp] 

After a SN explosion, the momentum of the ejecta is 
gradually transferred to the ISM during the expansion 
of the SN remnant until the final momentum-conserving 
snow-plough phase. As a result, SN-driving covers a wide 


^ This is a relatively rare event, as most SN remnants experience 
some interaction with dense gas before they reach a size of 30 
pc, and thus the compressive modes are usually not dominant, as 
explained in §5.1. 


The compressive power may exceed the solenoidal one at snap¬ 
shots with a SN remnant that expands in a uniform density region, 
making the baroclinic term negligible as e^lained below. This is 
the case of the specific snapshot of Figure [u 
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Fig. 9. — Left: Power spectra of solenoidal (solid lines) and compressive (dashed line) velocity components within a (100 pc)^ volume 
approximately centered around a SN explosion. The spectra are computed at three different times, to = 46.625 Myr, ti = 46.750 Myr and 
t 2 = 46.875 Myr. The time to is just before the SN explosion, while the times ti and t 2 are after the explosions, when the diameter of the 
remnant is approximately 40 and 80 pc, respectively. Between the times ti and t 2 , the expansion of the SN remnant brings the power in 
both the compressive and solenoidal modes to larger scales, as explained in the text. Right: Squared root of projected temperature (upper 
row of panels) and logarithm of projected density (lower row of panels) for the three times at which the power spectra are computed. At 
time ti the edge of the remnant is barely visible in the logarithm of projected density, while at time t 2 the gas density within the remnant 
has significantly decreased in the lowest density regions. 


range of scales and persists for a relatively long time 
after the explosion. This complex driving process can¬ 
not be primarily compressive because vorticity is read¬ 
ily generated by the baroclinic effect and amplified by 
the non-linear advection, as explained below. Even the 
ISM forcing immediately after the explosion is far from 
purely compressive, as hydrodynamical instabilities of 
the blast wave start already in the interior of the star, 
generating vorticity and making the ejecta very clumpy 
and asymmetr i c (e.g. Chevalier fc Klein||1978 Herant fc 
Woosley||1994[ IKifonidis et al.||2U06| |Cc)Uch et 


Wongwathanarat et al. 2015). In oiur simulation, the 
forcing at the moment of the thermal energy injection 
is not purely compressive either, because the accelera¬ 
tion from the corresponding pressure force has a non-zero 
curl (a baroclinic term) due to the non-uniformity of the 
medium (see the argument below). 

To examine the ratio of compressive and solenoidal 
modes, it is helpful to write down the equations of the 
flow divergence and vorticity. 


+ {diU^){djUi) = (Vp- Vp)/y - {Vh)/P, (1) 

and, 

Duj o 

— (cj • V)u + u;(V • u) = (Vp X Vp)/p . (2) 

The terms on the right-hand sides arise from the pressure 
term in the Navier-Stokes equation, and (Vp x 
in Equation (2) is called the baroclinic term, which con¬ 
tributes to generate vortical motions when the gradients 
of p and p are not aligned. The {lv • V)n term is known 
as vortex stretching and can amplify the vorticity. 

In our simulation, the flow velocity is driven by the 
pressure term. If we denote as Ps the contribution to 
the pressure from the SN thermal energy source, the ef¬ 
fective driving acceleration is — (Vps)/p- It follows im¬ 
mediately from Equations (1) and (2) that the diver¬ 
gence and curl of the effective acceleration are given by 
(Vp • Vps)/X - (V^Ps)/p and (Vp x Vps)/X, respec¬ 


tively. Because the SN locations are selected randomly 
in our simulation and due to the density fluctuations in 
the ambient gas, the pressure and density gradients (Vp 
and Vps) around the boundary of the initial SN sphere 
are not aligned in general, meaning that, at the instant 
of the SN energy injection, the effective driving acceler¬ 
ation at the boundary of the SN sphere is not curl-free. 
Therefore, the effective driving in our simulation is not 
purely compressive. 

If (Vp • Vps)/p^ dominates over (V^Ps)/p (which is 
supported by the fact that |Vln(p)| • |Vp| > |V^p| in a 
significant fraction of our simulated flow), the ratio of 
compressive to solenoidal power generated by the pres¬ 
sure source is determined mainly by the angle, 6>, be¬ 
tween the directions of Vp and Vps- Because in our 
simulation the SN locations are random, one may expect 
that the angle between the gradients is also random, so 
that ((cos6>)^) 1/3 and ((sin6>)^) 2/3. In that case, 

the effective driving generates the vorticity variance (cor¬ 
responding to solenoidal motions) twice faster than the 
divergence variance. Although of a qualitative nature, 
the above argument provides an explanation of why the 
solenoidal spectrum is typically larger than the compres¬ 
sive one. 

By monitoring solenoidal and compressive spectra as 
a function of time, we can observe that the expansions 
of the SN remnants bring both the compressive and 
solenoidal modes to larger scales, as shown in Eigure 
The power in compressive modes moves to larger 
scales according to the simplified ID model given earlier, 
which also explains the wavy features in the spectrum; 
the solenoidal power is transferred to larger scales by the 
expansion through the non-linear term in Equation (2), 
u;(V • u). Eigure^ shows the velocity power spectra in a 
(100 pc)^ volume centered around a SN explosion. The 
spectra are computed at three different times, one just 
before the SN explosion, and two after the SN remnant 
has reached an approximate size of 40 and 80 pc. The 
velocity field is multiplied by a tapered-cosine window 
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function before performing the Helmholtz decomposition 
in Fourier space, to gradually set the velocity to zero at 
the volume boundaries, making the velocity field peri¬ 
odic. The transfer of both solenoidal and compressive 
power to larger scale is clearly seen between the times ti 
and ^ 2 , as the remnant expands from 40 to 80 pc. 

At times or in regions not too close to very young 
SN remnants, the solenoidal and compressive modes 
have a chance to develop cascades toward small scales 
and to interact with each other, evolving toward a dy¬ 
namically and statistically relaxed state. The time- 
averaged spectra in Figure indicate that in the re¬ 
laxed state the solenoidal component of the velocity is 
larger than the compressive component at all k. The 
ratio, x{k) = Ec{k)/E^(k), is 0.2 at /c = 10. At 
the smallest wave numbers, the modes are almost in 
equipartition, x(/c) ^ 1. However, E^{k) is remark¬ 
ably shallow, while Ec{k) oc so the compressive- 

to-solenoidal ratio decreases rapidly towards larger wave 
numbers, x(/c) oc in the inertial range. Although 

they did not obtain clear power-law scaling and did not 
ident ify inertial-range slo pes in their SN-driven simula¬ 
tion, Balsara et al. (2004) also found that x(/c) decreases 
with increasing k and xW ^ 1 at large wave numbers. 

The different inertial-range slopes of Ec{k) and Es{k) 
and the rapid decrease of x(^) with increasing k are 
in contrast to previous studies, which typically found 
that x(^) is more or less constant in the inertial range. 
For example, simulations of weakly compressible turbu¬ 
lence with Mach numbers Atg ^ 1 showed that both 
compressive and solenoidal modes have Kolmogorov-like 
inertial-range sp ectrum, Edk) oc E^(k) oc and 

x{k) ^ 0.05 (e.g. Porter et al.||1998 , 1999, 2002). At the 
opposite extreme, simulations with purely compressive 
driving find x(^) ^ 1 in the inertial range, and the same 
Burgers-like slop es for both modes, Ec{k) ^ Es{k) oc k~‘^ 
(|Federrath|2013|) . Simulations with solenoidal (or mixed) 
driving and large Mach number. Als ^ 1. yield value s 
x{k) ^ 0.3 — 0.5 (Kritsuk et al. 2010 Federrath| 2013), 
with only a slight decline towards large wave nuuibers. 

An important difference between these studies and our 
simulation is that they all adopted a barotropic equa¬ 
tion of state, assuming either adiabaticity or isothermal- 
ity. The baroclinic effect is thus absent in all the sim¬ 
ulations mentioned above. As discussed earlier, when 
SNe explode in our simulation, the baroclinic effect from 
the pressure souree immediately drive solenoidal motions 
around the SN spheres at a rate similar to compress¬ 
ible modes. As the SN remnant and the flow evolve, the 
general baroclinic effect in Equation (2) is also likely to 
play an important role as the flow develops a cascade and 
evolves toward relaxation. As compressive motions in the 
form of shocks or expansions encounter a dense region 
in the flow, the baroclinic effect gives rise to shear and 
vortices around and within the region. Since the baro- 
clinic term depends on the gradients of the density and 
pressure, the conversion from compressive to solenoidal 
motions is expected to be more efficient at smaller scales. 
This explains why the compressive-to-solenoidal ratio de¬ 
creases rapidly with increasing k in our simulation. We 
thus argue that the existence of the baroclinic effect in 
our simulation is responsible for the different behaviors of 
Ec{k), Es{k) and x(^) compared to previous barotropic 


simulations. 

There has been recent interest in the regime of super¬ 
sonic turbulence with purely compressive driving, show¬ 
ing t hat it affects the probability distribution of gas den¬ 


sity ([Federrath et al.|2QQ8[ [Molina et al. ||2Q12|), the mag¬ 
netic held amplihcation by turbulence (4ederr ath et al. 


201 lat, and the inertial-range velocity scaling ([bchmidt 


et al. |2008 2009 iFederrathl 2013), with possible conse¬ 
quences tor models of star formation o r for turbulence 
in galaxy clusters (Porter et al. 2015). However, our 
result that SN-driving is not primarily compressive, par¬ 
ticularly at MC scales, suggests that features specific to 
isothermal flows with highly compressive driving may not 
apply to ISM turbulence. If SN shock-waves are the main 
energy source for MC turbulence, the driving accelera¬ 
tion would consist of a significant fraction of solenoidal 
modes, which arise through the baroclinic effect, when 
the SN shock impacts a cloud. Thus, idealized isother¬ 
mal simulations with purely solenoidal driving may bet¬ 
ter capture MC turbulence than isothermal simulations 
with purely compressive driving, and previous results on 
turbulent fragmentation based on solenoidal driving may 
not require significant correction. A careful study of the 
full implication for star formation of our result would 
require the evaluation of the compressive-to-solenoidal 
ratio specifically for the clouds formed in the SN-d riven 
turbulence, which we pursue in a separate work (|Pan 


5.2. Energy Injeetion Seale 

Using the velocity power spectrum, E(k)^ we can de¬ 
fine a length, Lin, that corresponds approximately to the 
scale where most of kinetic energy is contained, and can 
thus be interpreted as a characteristic scale of energy 
injection by SN explosions: 


_ 27r f k ^E(k)dk 
^ jE{k)dk ' 


(3) 


The time dependence of Lin is shown in Figure where 
we have also plotted the rms velocity, TEe value 
of Lin oscillates between approximately 50 and 100 pc, 
with many of the peaks corresponding also to peaks in 
ay. The time average and standard deviations are 70.5 
pc and 12.0 pc respectively. 

Figure also shows the transverse integral scale, L 22 - 
The longitudinal and transverse integral scales, Ln and 
L 22 , are defined as the integrals of the two-point velocity 
correlation functions: 


Ln = 




L22 = 




(4) 


where Rl and Rn are the longitudinal and transverse 
correlation functions. 


Rl{^) = {Ui{x)Ui{x^i)), 
Rn{E) — (pijiix^Ujiix -|- '^)), 


( 5 ) 


with U£ and Un being the velocity component parallel to I 
and one (of the two) transverse component perpendicular 
to respectively. Ln and L 22 are typically smaller than 
the injection length, Lin- In isotropic turbulence, exact 
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relations exist between Ln, L 22 and Lin- For exa mple, in 
incompressible turbulence, Ln = 2 L 22 = 3I/in/8 (Monin 


fc Iaglom]|1975 ). Our simulated flow is highly compress- 


ible, and we find that the longitudinal and transverse in¬ 
tegral scales are close to each other with (Fn) = 19.9 pc, 
and (1/22) = 19.4 pc. The near equality of Ln and L 22 
suggests that kinetic energies contained in solenoidal and 
compressive modes are comparable. Figure shows that 
at the scale of Lin the solenoidal and compressive spectra 
are in equipartition, in the sense that the solenoidal spec¬ 
trum is twice larger than the compressive one due to the 
extra degree of freedom. One can demonstrate that in 
such case Ln and L 22 should be exactly equal, as found 
in our simulation. Furthermore, in that case Ln and L 22 
are expected to be equal to Lin/4, which further explains 
the ratios, Lin/Ln = 3.5 and Lin/L 22 = 3.6. 

By comparing fo ur galactic-fountain simulations with 
different SN rates, |Joung et al.| (2009|) found that the 
energy injection scale decreases with increasing SN rate. 
In the model with approximately the same SN rate as in 
our simulation they obtained Lin = 87 pc. This value 
is approximately consistent with the one derive d here, if 
we account for the fact that Joung et al. (2009) assumed 
that 60% of the SN explosions are spatially correlated, 
as discussed above in §2.1, which enhances the forma¬ 
tion of super-bubbles and so should tend to increase the 
correlation length. _ 

Using their ga l actic- fountain simulation, |de Avill^ 
|& Breitschwerdt| (|2007|) computed the longitudinal and 
transverse integral scales, Ln and L 22 - They found that 
L 22 /L 11 = 0.5 — 0.6, consistent with the ratio in incom¬ 
pressible turbulence, implying that the overall compress¬ 
ibility of their simulated flow was likely low. Their mea¬ 
sured value of 75.2 pc for Lii would indicate an injection 
length scale of Lin — 200 pc. This injection length scale 
is significantly larger than Lin — 70 pc found in our sim¬ 
ulation. Based o n the dependence of the integral scale on 
the SN rate from Joung et al. (2009), this difference may 
be attributed to the much smaller SN rate in the simu¬ 


lations analyzed by de Avillez & Breitschwerdt (2007) 


6. VELOCITY STRUCTURE 

In order to characterize the turbulence, we compute 
the velocity structure functions, first for the whole com¬ 
putational volume, and then for individual clouds. The 
velocity structure functions of order p are defined as: 


Sp{£) = {\u{x + €) — u{x)\'^) oc 


(6) 


where the velocity component u is parallel (longitudinal 
structure function) or perpendicular (transversal struc¬ 
ture function) to the vector £ and the spatial average is 
o ver all pos it ions x . 

Boldyre\^ ( |2QQ2 ) proposed an extension to super- 


sonic tur bulence of the intermittency model b y She and 


Leveque ( [She fc Leveque|p~994| Dnbrulle||1994 ): 


C(p)/C(3)=p/9 + 1-(1/3UL 


(7) 


This velocity scaling has been found to provide a very ac¬ 
curate prediction for numerical simulatio ns of highly su¬ 
personic and super-Alfvenic turbulence ([Boldyrev et al. 


2QQ2| iPadoan et ^|2 QQ4b| |Pan & Scannapieco||2Ull|) 
Fadoan et al. showed that, as the rms IVlach 


u 

CL 



Fig. 10. — Time evolution of the energy-injection scale, Lin (thick 
solid line), computed from the integral of the inverse wave number, 
\Jk^ weighted by the velocity power spectrum, Eik) (see equation 
1^). The time-averaged value of 70.5 pc is shown by the horizontal 
dashed line. The dotted line shows the rms velocity, ctv (see the 
dashed line in Figure |m for the actual values of ctv). Almost all 
peaks in ctv correspona to peaks in L. The thin solid line shows 
the time evolution of the transverse integral scale, L 22 (see text). 
The average value, {L 22 } = 19.4 pc, is shown by the horizontal 
dashed-dotted line. 

number of the turbulence increases, the structure func¬ 
tion scaling varies from the She-Leveque scaling of in¬ 
compressible turbulence to the Boldyrev’s scaling, which 
has been interpreted as a gradual change of the Haus- 
dorff dimension of the most dissipative structures from 
1 (dissipation in filaments) to 2 (dissipation in sheets). 
Although the scaling of equation (|^ ma y not apply to 
flows driven by purely compressive forces (Schmidt et al.| 
2009), we have shown in §5.1 that the compressive modes 


are not dominant in SN-driven turbulence. 

Because of the limited extent of the turbulence inertial 
range in numerical simulations, the structure functions 
are usually power laws only over a very limited range 
of scales, if at all. Thus, the scaling exponents are usu¬ 
ally derived by normalizing the structure functions to 
the third-order one, which yields power laws that extend 
well into the (numerical or physical) dissipation range 
of scales. Thi s useful property is known as ‘extended 
self-similarity’ ( |Benzi et al.]|1993| ). 

The third-order structure function always yields a 
slope larger than unity in supersonic turbulence, while 
C(3) = 1 is an exact result in incompre ssible turbu¬ 
lence, the so called ‘4/5 law’ first derived by Kolmogorov 
(|1941|). This is because in supersonic turbulence the 
third-order structure functions should be computed with 
some density-weighting factor. For example, a density 
weight inspired by the assumption of constant energy 
trans fer of a compres sible flow, u r\j (£IpV^^ dFleckl 
1996), w as proposed by Kritsuk et al. (2QQ 7|) and further 
tested by Kowal & Lazarian| (|2UU7|) and byjbchmidt et al. 

(2008). In this work, we cornpute the structure functions 
either using the tracer-particle positions and velocities, 
or using the gas velocity values on a uniform mesh, with 
a density-weighting method that is equivalent to the use 
of tracer particles. This densi ty-weighting method is dif¬ 
ferent from that proposed by Kritsuk et al. (2007). 

To calculate the structure functions based on tracer 
particles in a MC, we search particle pairs at given dis- 
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Fig. 11.— Density-weighted longitudinal velocity structure func¬ 
tions, (see equation (|^), of orders p = 1, 2 and 3, plotted 

versus the separation i, computed over the whole computational 
volume (filled circles). The plots show the time average of the 
structure functions from 24 snapshots covering uniformly the time 
interval t = 33 — 56 Myr (one snapshot per Myr). The solid lines 
corresp onds to the structure function slopes predicted by |Boldyrev| 
(|2002|) and con firmed in simulations of randomly-driven i sother¬ 
mal turbulence (iB oldvrev et al.|2002[ iPadoan et al.|2004b| IParT&l 
|Scannapieco|201l|), (,(,-*-) = = U.Y4, (,(3) = l.U. Allliougli 

TTTe^r^rovlHelirrexcellent fit, these solid lines are not obtained by 
fitting the data. 

tances, compute their relative velocities, and average 
over all the particle pairs to obtain the relative velocity 
moments at different orders, p. For example, the p-th 
order structure function is computed as 

1 ^ 

(8) 

n=l 

where N is the number of pairs separated by a distance 
of and Uin — U 2 n is the relative velocity of the n-th 
pair of particles. These tracer-based structure functions 
have a built-in density weighting, because the number 
of particle pairs at a given distance depends on the flow 
densities at the particle positions (the tracer particles are 
initialized with a number density proportional to the lo¬ 
cal gas density). For structure functions over the entire 
simulation box, we compute grid-based structure func¬ 
tions with a density weighting defined as 

cdw/M _ {p{x)p{x + l)\u{x + t)- u{x)\P) 

{p{x)p{x + e)) • 

It is straightforward to prove that Sp{£) and are 

statistically equivalent to each other, if the number of 
tracers is large enough to sufficiently sample the flow 
density field. To show this, consider two infinitesimal 
volumes, dVi and dV 2 , around two points, 1 and 2 , sep¬ 
arated by a distance of i. When computing the tracer- 
based structure functions, these two volumes would con¬ 
tribute to a total number of N 12 = {fip/p)‘^pip 2 dVidV 2 
particle pairs at a distance of where rip and p are the 
mean particle number density and mean flow density, re¬ 
spectively. In other words, the two points 1 and 2 are 
essentially counted W 2 pip 2 times in the computa¬ 
tion of This is equivalent to a density-weighting 

factor of (X pip2^ suggesting that S^p{t) is equal to the 


Fig. 12. — Same as in Figure El but for a single snapshot at 
t = 44.5 Myr. 

grid-based structure function, 

The density-weighting scheme adopted here is of par- 



{p{x)p{x -I- i)Ui{x)Ui{x -I- i)Uj{x -I- 

{p{x)Uj{x)p{x -^i)) oc ij, (10) 

where p is the pressure of the flow. In highly supersonic 
turbulence, the pressure term may be neglected, and we 
have 

{p{x)p{x i)Ui{x)Ui{x i)Uj{x i)) (xij. (11) 

The quantity {p{x)p{x + i)ui{x)ui{x -f- i)uj{x i)) is 
closely related (although not exactly equivalent) to the 
density-weighted third order structure function, £' 3 ^. We 
therefore expect a linear scaling for our density-weighted 
third order structure function, S' 3 ^ oc which turns out 
to be confirmed by our simulation for both the entire flow 
and individual MCs. 

6 . 1 . Global Velocity Scaling 

We first compute the velocity structure functions aver¬ 
aged over the whole computational volume. The velocity 
field is first remapped onto a uniform mesh of 512^ cells, 
and then the structure functions are computed with the 
density-weighting method described above. The volume 
contains large voids of very low density, with a very low 
number of tracer particles, so we prefer this mesh-based 
method instead of the equivalent tracer-particle method 
(which we use below for the structure functions inside 
MCs). To speed up the calculations, we only consider 
velocity differences along the three main orthogonal di¬ 
rections, and 16 values of cell distances, £. Even with 
this limitation, the sample size is large enough to yield 
reliable low-order statistics. 

Figure pd] shows the first, second, and third order longi¬ 
tudinal velocity structure functions, 5'p^(-^), plotted ver¬ 
sus the separation computed over the whole volume 
as described above, and time-averaged using 24 snap¬ 
shots covering uniformly the time interval t = 33 — 56 
Myr (one snapshot per Myr). They are well approxi¬ 
mated by power laws in the approximate range of scales 
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Fig. 13.— Structure function slopes versus time, for structure 
functions computed in each of the 188 snapshots in the time interval 
t = 33 — 56 Myr. The horizontal dotted lines show the values of the 
first or der slope, C(l). predic ted by the She-Leveque intermittency 
model (iShe fc L eveque|1994|> (lower line) and the Boldyrev model 
(|Boldyrev|2UU2| (upper linej. The horizontal solid lines correspond 
to tlie predictions for ^(2) by the same models, and the horizontal 
dashed line corresponds to C( 3) = 1, an exact ma thematical result 
for incompressible turbulence |Kolmogorov| ( |l941| >. 

between 3 and 30 pc. The values of the exponents are 

c(l) = 0.39T0.01, C( 2 ) = 0.75T0.01, C(3) = 1.13T0.02. 

To compare with values previously found in studies 
of randomly-driven supersonic isothermal turbulence, in 
this and following figures we over-plot lines showing the 
slopes from eq. Notice that eq. [^only gives the expo¬ 
nents of the velocity structure functions without density 
weighting and normalized to the third-order one. Thus, 
in this comparison, we make the reasonable assumption 
that the exponents normalized to the third order should 
be the same for the unweighted structure functions as for 
the density-weighted ones. However, future works should 
recompute the velocity structure functions of randomly- 
driven, supersonic isothermal turbulence simulations us¬ 
ing the same density-weighting scheme as in this work. 

While the first and second order structure functions in 
Figure have exponents within 3 sigma and 1 sigma 
respectively from the values of equation 0 C(3) is sig- 
nificantly larger than unity, despite the time averaging. 
We find that this deviation from unity, and even from a 
power law, occurs in snapshots during periods with the 
highest SN rate, when it is more likely that a snapshot is 
very close in time to a very recent SN explosion. During 
the brief, initial period of the SN bubble expansion, SN 
driving has a direct effect on a range of scales. On the 
contrary, in snapshots that are not too close in time to 
SN explosions, the SN remnants have already expanded 
to large scale, the turbulence has had time to relax, and 
the third order structure function is found to be a nearly 
perfect power law with C(3) = 1. Figure shows an ex¬ 
ample of velocity structure functions from a single snap¬ 
shot, at t = 44.5 Myr, that is not too close in time to 
a SN explosion. The third order structure function is 
now a power law, with the expected numerical decay be¬ 
low approximately 3-4 pc, and the third order slope is 
indistinguishable from unity. The exponents in this sin¬ 
gle snapshot are C(l) = 0.41 ± 0.01, C(2) = 0.77 ± 0.02, 
C(3) = 1.00 ±0.02, consistent with equation [^within one 
or two sigma. 


Fig. 14. — Velocity structure functions obtained from the average 
of the structure functions of the 15 most massive clouds selected 
at time t = 34.2 Myr. The structure functions of the clouds are 
computed from the posit ion and vel ocity of their tracer particles 
(see text). As in Figures [FI and |12[ the solid lines are the model 
predictions for the structure function slopes, not fits to the data. 


To better quantify the uncertainty of these exponents, 
and to illustrate their time dependence, we compute their 
values for each snapshot in the time interval t = 33 — 56 
Myr (8 snapshots per Myr), by fitting the structure func¬ 
tions as a function of i in the range between 3 and 
30 pc, and plot them versus time in Figure The 
two pairs of solid and dotted horizontal lines sEow the 
predicted values for incompressible (lower values) and 
supersonic (higher values) turbulence (She-Leveque and 
Bodyrev scaling respectively). One can see that the val¬ 
ues are usually within the predicted ones for the first and 
second order exponents, while the third order exponent 
is systematically above unity during the first and last 
thirds of that time interval, when the SN rate and the 
kinetic energy appear to be a bit higher than in the mid¬ 
dle third (see Figured. The time average of the single¬ 
snapshot values are C(T) = 0.39±0.04, C(2) = 0.75±0.05, 
C(3) = 1.1 ±0.2, well within one sigma of the values from 
equation 

The logarithmic fluctuations in Figure (the y-axis 
of the plot is logarithmic) are significantlyTarger for the 
third order than for the first and second orders, mean¬ 
ing that deviations from the expected scaling due to the 
effect of SN driving is stronger for higher order statis¬ 
tics. Therefore, it is not convenient to take advantage 
of extended self-similarity and measure the exponents by 
plotting the structure functions versus the third order 
one, instead of versus at least not for the first and sec¬ 
ond order case, as their time evolution would then have 
much larger fluctuations, due to the larger fluctuations 
of the third order one. 


Jonng fc Mac Low (2006) and de Avillez & Breitschw- 
erdt| ( |2007[ ) computed the velocity scaling exponents from 


galactic-fountain simulations and found a scaling law 
consistent with Boldyrev’s prediction. However, because 
of insufficient dynamic range below the driving scale, 
their simulations did not yield a power-law inertial range, 
and so the exponents could be computed only relative to 
the third order one (taking advantage of the extended 
self-similarity). Therefore, their result cannot be com¬ 
pared directly with the ISM velocity scaling derived from 
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Fig. 15. — Same as in Figurebut for a single MC, the seventh 
most massive one among those selected at time t = 45.6 Myr (MC 
6 in the top panel of Figure]^. 

observations. Furthermore, their simulations did not 
reach the necessary spatial resolution to describe the evo¬ 
lution of individual SN remnants and their interaction 
with MCs in detail, and to study the velocity scaling 
within individual clouds, to test if MC turbulence is con¬ 
sistent with SN driving. 

Because our simulation yields power-law velocity 
structure functions as a function of and the values 
of the exponents are the same as in previous numerical 
studies of MC turbulence based on random driving with 
a large-scale volume force, we conclude that the use of 
an artificial force in those previous studies did not result 
in incorrect velocity scaling, so no major corrections to 
IMF and SFR models based on turbulent fragmentation 
should be needed. However, because deviations from the 
average scaling laws are found in snapshots that are very 
close in time to SN explosions, it may be worthwhile 
investigating the possibility of minor effects on star 
formation resulting from such deviations. 


6.2. Velocity Scaling within MCs 

In order to test if SN driving can generate the same 
velocity scaling also inside MCs, despite their density 
contrast, we compute the velocity structure functions in¬ 
side MCs selected from our simulation as described in §3. 
To better constrain the scaling exponents, we have com¬ 
puted the velocity structure functions for the 15 most 
massive clouds of each snapshot. Because of the very 
complex cloud shapes, we find it convenient to compute 
the structure functions based on the position and veloc¬ 
ity of the tracer particles, Sp^{i) (see eq. MCs are 
regions of density enhancement, so their velocity field is 
sampled well by the tracer particles. In fact, they contain 
so many tracer particles that, to speed up the calculation, 
we randomly select a number of particle pairs 500 times 
smaller than all possible pairs in each cloud, resulting 
in approximately 2 to 200 million pairs per cloud. As a 
further simplification, we also compute the differences of 
each velocity component irrespective of their orientation 
relative to the separation vector, so we do not dis¬ 
tinguish between transversal and longitudinal structure 
functions. 

We find velocity differences growing with scale up to 


Fig. 16. — Same as in Figurebut for the fifth most massive 
MC among those selected at time t = 56.1 Myr (MC 4 in the bot¬ 
tom panel of Figure!^. This is an example of structure functions 
clearly affected by deviations at small scales due to the effect of 
nearby and recent SN explosions. During the brief period of time 
of the early expansion of a SN bubble, excess energy is found at 
small or intermediate scales relative to the ‘undisturbed’ structure 
functions. 


10-30 pc, depending on the cloud size. As an example. 
Figure shows the average of the structure functions 
of the 15 most massive clouds selected at time t = 34.2 
Myr. The structure functions inside these GMCs are con¬ 
sistent with the global ones and are well approximated by 
power laws down to a scale of approximately lOdx = 2.4 
pc. At a separation £ = Adx = 0.96 pc, the velocity 
differences are only approximately 30% below the value 
extrapolated from the inertial-range scaling, as shown by 
the first order structure function in Figure 

In Figures and we show examples of structure 
functions of individuaTMCs, one with no evidence of the 
direct effect of SN driving (Figure [T^, and one with sig¬ 
nificant deviations at ^ < 1 pc, due to a recent nearby 
SN explosion (Figure [l^. 

These results show that, despite the large contrast be¬ 
tween the MC density and the average ambient density, 
the kinetic energy of SN explosions is effectively trans¬ 
ferred into turbulence within individual clouds, where 
it establishes the usual scaling laws of supersonic turbu¬ 
lence, all the way to the smallest scales where the simula¬ 
tion is affected by numerical dissipation. Real MCs also 


exhibit power-law velocit y scaling (e.g. Heyer & Brunt 
2004 Padoan et al. pop'd ). The slo pe of C(^) = 0. 8 =b 0.1 


derived for the Pers eus region by 


Padoan et al. (2006) 


using the method by Lazarian & Pogosyan (|2000 ) is tor- 
mally consistent with the scaling laws derived here for the 
MCs of our simulation. On th e other hand, the princ iple- 
component-analysis results by Heyer & Brunt (2004) give 
a very large slope, ({2) = 1.12 ± 0.04, which is hard to 
interpret, as it is steeper than the scaling from the Burg¬ 
ers equation that models an infinitely compressible flow. 
A careful study of the consistency between the structure 
function slopes of the SN-driven simulation and of the 
observations requires the computation of synthetic ob¬ 
servations and is beyond the scope of this work. 

All the plots in this section adopt an arbitrary 
normalization of the structure functions. The actual 
normalization of MC turbulence, that is the velocity 
dispersion of clouds of a given size, is discussed below. 
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Fig. 17.— Virial parameter versus energy ratio for clouds selected 
from three snapshots of the simulation prior to the inclusion of self¬ 
gravity and covering a time span of 4 Myr. The solid line shows 
the average ratio of 1.20 between avir and 2E]^lEg. The dotted 
line marks the maximum value of that ratio. 

comparing the velocity-size Larson relation of our clouds 
with that from the observations. 


7. VIRIAL PARAMETER AND CLOUD STRUCTURE 



2EJE^ 


Fig. 18. — Virial parameter versus energy ratio for clouds selected 
from three snapshots after the inclusion of self-gravity and covering 
the last 4 Myr of the sim ulation. The solid and dotted lines are 
the same as in Figure pT| showing that the average ratio of avir 
and 2Ek/Eg and its scatter are nearly unchanged for most clouds, 
relative to their value before the inclusion of self-gravity. Clouds 
with very low values of 2E'k/Eg contain collapsing cores whose 
gravitational energy has been included in the computation of the 
cloud Eg. They are star-forming clouds with relatively low value 
of avir, but do not show evidence of global collapse. 


Larsonj (1981) interpreted the MC velocity-size rela- 
tion he discovered as due to a turbulent cascade in the 
ISM. Because of the very large Reynolds number of the 
observed motions in the cold ISM, the development of 
a turbulent cascade is unavoidable, and both analytical 
arguments and numerical simulations have demonstrated 
that Larson relations can be viewed as the natural conse - 
quence of supersonic turbulence (e.g. Kritsuk et al.|2013 ). 
Nevertheless, the velocity-size relation has also been in¬ 
terpreted as the consequence of the MC se lf-gravity (e.g. 


Solomon et al. 1987 Heyer et al. 2009), because MC 
virial masses are often comparable to tne masses esti¬ 
mated from the CO luminosity. 

The relative importance of turbulence and self-gravity 
is measured by the vir ial parameter, introduced by 


Bertoldi & McKee (1992): 


_ ^ 40 / iff y 2Ek 

“ GMei 37r2 Ldyn; ’ 


( 12 ) 


where is the one-dimensional velocity dispersion, and 
the dynamical time is defined as: 


^dyn — -Rci/<^v,3D* 

The last equality in ( [T^ is exact in the case of an ideal¬ 
ized spherical cloud oiuniform density. For more realistic 
cloud mass distributions, the virial parameter is only an 
approximation of the ratio of kinetic and gravitational 
energies. 

To estimate the relative importance of turbulence and 
self-gravity in clouds from our simulation, we have com¬ 
puted the virial parameter and the kinetic and gravi¬ 
tational energies of clouds selected from six snapshots, 
three before and three after the inclusion of gravity in 
the simulation, using a threshold density nH,min = 100 
cm“^. The virial parameter of a cloud is measured using 
the positions and velocities of the tracer particles belong¬ 
ing to that cloud, and defining the cloud radius as the 


rms of the particle positions: 

3 N 


Rc\ = 




i=l n=l 


1/2 


(14) 


where Xi = ^i,n/N are the components of the 

mean particle position (the cloud’s barycenter) and N 
is the total number of tracer particles in the cloud. 

While the virial parameter only depends on global MC 
properties (mass, radius and rms velocity), the ratio of 
kinetic to gravitational energy is sensitive to the cloud 
internal structure (mass distribution and sha pe) and to 
correlations bet ween density and velocity (e.g. Federrath 
fc Klessen|2Q12 ). Thus, the comparison between and 
2F/’k/L/g is a useful tool to probe the internal structure 
of MCs and its evolution under the effect of self-gravity. 
As in the case of the virial parameter, we compute E\^ 
and Eg of a cloud using the velocities and positions of 
the tracer particles in that cloud: 


II 

to 1 

A 

(15) 

N N . 

^ in 

(16) 


i=i j=i 

where N is the total number of tracer particles in the 
cloud, m is the mass associated to a tracer particle, Ui is 
the modulus of the velocity of the i-th particle, and rij is 
the distance between the i-th and the j-th particles. This 
expression for Eg assumes that the cloud is isolated, as 
in the definition of the virial parameter. In future work, 
this expression should be contrasted with a formulation 
that accounts self-consistently for the surrounding mass 
distribution by using th e actual gravitational poten tial 
from the simulation (e.g. Federrath fc Klessen||2Q12 ). 
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Fig. 19.— Probability distributions of the viral p aram eter for the 
same three snapshots before self-gravity as in Figure |17| (shaded his- 

a ram) and for the three snapshots after self-gravity as in Figure 
(unshaded histogram). 


Fig. 20.— Probability distributions of mean clou d de nsity for the 
same three snapshots before self-gravity as in Fignre |17| (shaded his- 

S ram) and for the three snapshots after self-gravity as in Figure 
(unshaded histogram). 


Figure shows the comparison between and 
2E\^lEg for clouds selected from three snapshots cover¬ 
ing a time span of 4 Myr before the inclusion of self¬ 
gravity. Figure 18 shows the same plot, but based on 
three snapshots after the inclusion of self-gravity, cover¬ 
ing the last 4 Myr of the simulation. On average, the 
clouds of Figure [l^ are selected 9 Myr after the inclusion 
of self-gravity, while their average free-fall time, based 
on the density estimated as nH,ci = Mci/(47ri?y3) (see 
Figure 20), is 4.4 ± 2.1 Myr, and their average dynam¬ 
ical time, tdyn = 7?ci/crv,3D, is 2.6 ± 1.5 Myr. Thus, 
self-gravity has been active for approximately two cloud 
free-fall times and four cloud dynamical times on aver¬ 
age, with significant changes occurring only in regions 
that have too small filling factors (e.g. small collapsing 
cores) to change global statistics significantly. Indeed, 
gravity is not expected to be important for clouds with 
Q^vir > 1- Selecting clouds with a^ij. < 1.0, 0.5 and 0.25, 
we find a mean free-fall time of 2.1 ± 1.3 Myr, 1.5 ± 1.3 
Myr and 1.2 ± 1.0 Myr respectively. Thus, for the clouds 
where it should be most important, self-gravity has been 
active for 4 to 8 free-fall times on average. To estab¬ 
lish whether the absence of significant trends during this 
time interval will continue for the lifetime of the various 
structures requires, as discussed in Section 2.1, future 
simulations that include the supernova feedback from the 
massive stars produced by the simulation itself. 

Figure shows that a^ir provides a remarkably good 
approximation of the energy ratio, 2E\^lEg^ despite the 
complexity of the cloud structure. We find the ratio 
ayir/{2E\^/Eg) = 1.20 on average, constant over three 
orders of magnitude in E\^/Eg (we have verified that it 
is constant also over the full range of cloud masses) and 
with a small scatter (the standard deviation is approx¬ 
imately 20% of the mean); it is also nearly unchanged 


after gravity is included in the simulation. Figure 18 


shows that most clouds follow approximately the same 
relation of versus 2E\^lEg as in the case without 
self-gravity, except that a few of them have significantly 
lower values of 2E\^IE^^ because they contain collapsed 
cores that have also been included in the computation of 
the total E^. The comparison of the two figures, without 
and with gravity, shows the main effect of self-gravity 


is to cause the collapse of dense cores inside the clouds 
(hence star formation), while the cloud virial parameter 
is not strongly affected. 

In Figure we show the probability distribution of 
for the three snapshots without gravity (shaded his- 


a 


togram) and with gravity (unshaded histogram), where 
we have included also the star-forming clouds with very 
low values of 2E\^lEg. The two distributions are very 
similar, with the one including gravity slightly shifted to¬ 
wards lower values; the mean values are 8.5 without grav¬ 
ity and 6.6 with gravity. The small shift between the two 
distributions shows that self-gravity causes some amount 
of cloud contraction, but not a significant change in 
global cloud structure. This is further confirmed by the 
histograms of cloud density shown in Figure where 
the cloud density is defined as nH,ci = 

The histogram for the clouds with self-gravity is only 
slightly shifted to higher density, with the mean value 
changing from 183 cm“^ to 264 cm“^, before and after 
the introduction of gravity respectively. This small in¬ 
crease in cloud density shows that self-gravity does not 
cause a global cloud collapse, even if it drives star for¬ 
mation through the collapse of dense cores within MCs. 
Only clouds with ^ 10 contribute to star forma¬ 
tion, according to Figure with most star formation 
occurring in clouds with < 3, in agreement with re¬ 
cent studies of the SFR in supersonic turbulence, show¬ 
ing that the SFR is mainly controlled by the virial pa¬ 


rame t er (Padoan fc NordlundpQll Federrath & Klessen 
2012 Fadoan et al. '2012^ Future simulations, where 
self-gravity is active for much longer than the 11 Myr 
period of our run, and where selfconsistent supernova 
feedback is included, will be needed to further test the 
above results. 


8. MC LIFETIMES 

While the virial parameter estimates the relative im¬ 
portance of turbulence and self-gravity, the comparison 
of the cloud lifetime with either the dynamical time or 
the free-fall time provides a definitive assessment of the 
actual dynamical state of clouds. For example, short¬ 
lived clouds are evidently not gravitationally bound, 
while their virial parameter may be of order unity and 
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thus would not allow to draw a definitive conclusion 
about their dynamical state. 

Being so difficult to constrain observationally, MC life¬ 
times are a valuable outcome of numerical modeling. 
They can be measured in our simulation thanks to the in¬ 
troduction of tracer particles, but only up to a maximum 
age of 23 Myr, the time interval with tracer particles in 
our simulation (11 Myr if we considered only clouds dur¬ 
ing the time interval with self-gravity). Although this 
maximum age is comparable to or smaller than the life¬ 
time of the largest clouds, it is at least significantly longer 
than the mean free-fall time (4.4 Myr) and the mean dy¬ 
namical time (2.6 Myr) of the clouds in the simulation, 
so it allows us to evaluate the influence of gravity on the 
cloud lifetimes for most clouds in our study. 

Dobbs and Pringle (2013) measured MC lifetimes by 
defining the continuation of a cloud at a later time as 
that with the largest number of particles in common with 
that cloud. The time when the number of particles (and 
the mass) in common drops to half or less of that in the 
original cloud marks the end of the cloud lifetime. Cloud 
precursors and formation times are defined in the same 
way, by considering the clouds at earlier times with the 
largest number of particles in common with the original 
cloud. As pointed out by the authors, this method has 
the drawback that clouds precursors or continuations are 
often not found with more than half of the original mass 
because of changes in density contours, rather than a real 
cloud dispersion. For example, based on this method, a 
cloud may seem to have dispersed after a certain time, 
but may later reappear. 

After verifying in our own simulation that this lifetime 
definition is indeed very uncertain, we have chosen to 
estimate cloud lifetimes with a different method that is 
independent of the specific density contours of clouds in 
past and future snapshots. Because the formation of a 
cloud implies converging flows (even in the absence of 
self-gravity), and its dispersion requires diverging flows, 
we simply define the lifetime of a cloud as the time in¬ 
terval during which the cloud radius, defined always by 
all the tracer particles belonging to that cloud, is within 
a factor of two from the radius at the time the cloud is 
selected. In other words, we follow the cloud tracer par¬ 
ticles in the past, until their radius has doubled in size, 
which marks the formation time of the cloud, and in the 
future also until the radius has doubled in size, which 
marks the dispersion time of the cloud. 

Figure!^ shows the formation time, tform (interval be¬ 
tween cloud formation and time of cloud selection), and 
dispersion time, tdisp (interval between time of cloud se¬ 
lection and cloud dispersion) versus the cloud mass, using 
all the clouds more massive than approximately 700 Mq 
from each of 12 snapshots covering the whole time in¬ 
terval with tracer particles (the time between snapshots 
is 2 Myr). The sum of the two times gives the cloud 
lifetime, tufe = tform + ^disp, but we have plotted the 
two times separately (empty squares and circles) because 
in many cases we can identify the cloud formation time 
and not its dispersion time (for clouds selected towards 
the end of the simulation), or vice-versa (for clouds se¬ 
lected shortly after the introduction of tracer particles). 
The average values of tform and tdisp are biased by the 
large number of upper limits (not plotted to avoid con¬ 
fusion). In order to obtain a nearly unbiased estimate 



Fig. 21.— Dispersion time (empty circles) and formation time 
(empty squares) for clouds selected from 12 snapshots of our sim¬ 
ulation, uniformly distributed at intervals of 2 Myr, to cover the 
whole time interval with tracer particles. Filled circles mark the 
values of dispersion time of clouds selected from the first snapshot 
of the series; filled squares mark the values of formation time of 
clouds from the last snapshot. 


of the average values, we consider only clouds selected 
from the first and the last snapshots of the series (filled 
circles and filled squares, respectively, in Figure 21). The 
first snapshot has 12 clouds more massive than^O Mq, 
yielding 10 measured values of tdisp and only two upper 
limits; the last snapshots has 40 clouds more massive 
than 700 Mq, 39 of which with a measured value of tform 
and only one with an upper limit to tform- Based on 
these measurements alone, we find (tdisp) = H-l Myr 
and (tform) =9.2 Myr. By assuming that clouds are se¬ 
lected at a random moment of their lifetime, we should 
expect (tiife) = 2 (tform) = 2 (tdisp), and so we can aver¬ 
age together all the 49 values and obtain (tufe) = 21.4 
Myr. 

This lifetime is based on clouds covering over two or¬ 
ders of magnitude in mass, and only on 39 measurements. 
We can further improve our estimate of MC lifetimes by 
using all measurements and, at the same time, derive 
the mass dependence of the lifetimes, if we normalize 
the lifetime to the cloud dynamical time, tdyn- Figure 


22 plots tiife versus tdyn of 64 clouds for which we have 


measured values of both tform and tdisp (filled circles), 
besides the formation and dispersion times for all the 
other clouds where these are measured (er npty squares 
and circles respectively). The plot in Figure!^ shows an 
approximately linear correlation between tufe and 4 tdyn • 
The mean and standard deviation of the ratio of the two 
times are 

hifeAdyn = 4.5 ± 1.4. (17) 


As shown by Figure this estimate is derived from 
clouds with Mci < 6 x 10^ Mq, and thus it should be 
considered only as an extrapolation when applied to more 
massive MCs. Nevertheless, this result for the tiife/tdyn 
ratio is consistent with the measured formation times 
of the most massive MCs in the simulation. Figure 
shows that the four most massive GMCs, with masses 
around 10^ Mq, have (tform/^dyn) close to two, which 
would imply (tiife/tdyn) ~ 4, in the absence of selection 
biases due to the limited time interval covered by the 
tracer particles. 
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Fig. 22. — Cloud lifetime versus cloud dynamical time for clouds 
with measured dispersion and formation times (filled circles). 
Clouds with only dispersion or formation times are plotted as 
empty circles and empty squares respectively. The dashed line 
corresponds to the long-dashed line to = dt^yn- 


Furthermore, while we cannot use all the values of 
tform and tdisp from Figure to estimate an unbiased 
average lifetime from (tufe) ~ 2(tform) ~ 2(tdisp), we 
can still use all of the corresponding ratios, tform/^dyn 
and tdisp/^dyn, to estimate an unbiased average ra¬ 
tio of lifetime to dynamical time from (tiife/tdyn) ~ 
2(tform/tdyn) ~ 2 (tdisp/^dyn) , if wc assumc that this ra¬ 
tio is independent of cloud mass. Indeed, the dashed 
and dotted lines in Figure 23 show that (tform/tdyn) ~ 
(tdisp/tdyn) ~ 2, using valucs over the whole mass range, 
consistent with the estimate (tiife/tdyn) ~ 4, based on 


clouds with Mci ^ 6 x 10^ Mq. 

This result shows that both the formation and the dis¬ 
persion of the MCs in our sample take two dynamical 
times, on average. This is an indication that both the 
formation and the dispersion of the MCs in our sample 
is controlled by the turbulence, with little influence of 
self-gravity. Because of the non-negligible scatter in the 
ratio of cloud lifetime to dynamical time, one may expect 
that at least the clouds with the largest ratios may have 
longer lifetime due to their self-gravity. This is not the 
case: we have verified that there is actually a positive cor¬ 
relation between tiife/tdyn and avir,meaning that larger 
values of tiife/tdyn are usually due to smaller values of 
tdyn because of larger (hence larger avir), rather than 
longer tufe as a consequence of a lower Thus, there 
is no significant imprint of self-gravity in the cloud life¬ 
times, even if more than half of our clouds are selected at 
a time after self-gravity has been included in the simula¬ 
tion. Future simulations, where selfconsistent supernova 
feedback allows longer runs with selfgravity, are needed 
to test if this lack of significant imprint of selfgravity con¬ 
tinues, and if it extends to MCs with longer lifetimes and 
to higher surface density MCs. 

We should also stress the caveat that, for the most mas¬ 
sive MCs of ^ 10^ Mq we could only measure formation 
times, and not dispersion times, due to their long lifetime 
(and dynamical time) and the limited duration of the 
simulation. Thus, we cannot rule out that, at least the 
most massive MCs, could have dispersion times signifi¬ 
cantly longer than two dynamical times. However, that 
would imply dispersion times longer than 20 Myr (life¬ 


Fig. 23.— Ratio of lifetim e to dynamical time versus mass for 
the same clouds as in Figure The long-dashed line shows the 
mean ratio for the clouds with both dispersion and formation times 
measured (filled circles), {(tform + ^disp)/^dyn} = 4.5. The dotted 
and the dashed lines the mean ratio for clouds with only formation 
or dispersion times respectively, {tform/^dyn} = 2.2, (tdisp/^dyn) = 
2.4. All these average values are consistent with (tiife/^dyn) ~ 4, 
independent of cloud mass. 


times longer than 40 Myr) for such clouds, a timescale 
over which the extra energy injection from SN explo¬ 
sions of locally formed massive stars would presumably 
succeed in dispersing the clouds, even if the general ISM 
turbulence could not (see discussion in §2.1). 

To derive actual values of cloud lifetime as a function 
of cloud mass, taking advantage of our result (17), we 
can use the expression (25) for the average cloud ^nam- 
ical time derived in §l(kT^ which gives an average cloud 
lifetime of 


tute = 22.5 Myr (Md/lO^ (18) 


9. MAGNETIC FIELD IN MCS AND MC FORMATION 

Our simulation adopts a mean magnetic-field strength 
consistent with the Galactic one (see §2), so the magnetic 
field inside clouds selected from the simulation should be 
comparable to that in real MCs. We have already shown 
in Figures]^ andthat the mean magnetic energy is not 
far from equipartition with the mean thermal and kinetic 
energies averaged over the whole volume, while the en¬ 
ergy ratios are much larger in the dense gas. This clear 
energy separation in dense gas, with (^k,d/^m,d) = 25.1 
and (^m,d/^th,d) = 9.8, is the necessary consequence 
of the near equipartition at the largest scales. Being 
only mildly super-Alfvenic, large-scale compressive mo¬ 
tions cannot compress the mean magnetic field by a large 
factor, so the density enhancement of MCs is largely 
achieved with compressions along field lines, resulting in 
a mean magnetic field strength in the dense gas not much 
larger than the total mean field. The mean magnetic field 
of 4.6 jiiG is amplified by the SN-driven turbulence to an 
rms value of 7.2 /iC, averaged over the whole volume and 
between t = 33 and 56 Myr. The rms field strength in 
the dense gas is 12.8 /iC, not even a factor of two larger. 

To investigate the role of the magnetic field in individ¬ 
ual MCs, we consider the same catalog of 1,547 clouds as 
in the comparison with the observations discussed in the 
next section. The clouds are selected from 7 snapshots 
during the final 6 Myr of the simulation, at a resolution 
of 0.49 pc and with a density threshold of nH,min = 200 
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Fig. 24.— Magnetic field strength versus clou d mass for the same 
sample as in the Larson relations of Figures [m] and [M] The dashed 
line shows the mean magnetic field averaged over the whole com¬ 
putation volume (also the initial mean field). Empty circles corre¬ 
spond to the mean value of the magnetic field of all tracer particles 
in each cloud, while filled circles give the rms value. 

cm“^. We compute both the mean and the rms magnetic 
field of each cloud using the values sampled by the tracer 
particles, (B) = and , 

where Bi is the magnetic field strength sampled by the 
particle i in a given cloud, and N is the total number of 
particles in that cloud. These magnetic field values are 
plotted versus cloud mass in Figure where the hori¬ 
zontal dashed line represents the mean magnetic field in 
the computational volume, Bq = 4.6 /iG. The mean field 
in the clouds is approximately 10 /iG on the average, only 
twice larger than Bq^ and independent of cloud mass. We 
have verified that the mean magnetic field strength of the 
clouds is also independent of their mean gas density. 

The relatively small increase of the cloud mean mag¬ 
netic field relative to Bq and its independence of gas den¬ 
sity are characteristic of trans-Alfvenic su personic tur- 

and further 


1999) 


bulence (Padoan & Nordlund 1997 
illustrates that IVlCs must be formed by compressive mo¬ 
tions primarily along magnetic field lines, due to the 
non-negligible magnetic pressure prior to the compres¬ 
sion and cooling of the low-density gas. As the gas is be¬ 
ing compressed into a nascent MG by random large-scale 
motions, the increasing density and decreasing cooling 
time cause a drop in both the Alfven and sound speeds 
( Padoan et al.|[2CilQ ). As a result, the turbulence within 
a IVIC is super-Altvenic and highly supersonic, while 
the larger-scale flows responsible for its formation are 
trans-Alfvenic and mildly supersonic. Because in super- 
Alfvenic turbulence the magnetic field is amplified by 
compressions, as shown by a posit ive B — n correlation 
(Padoan & Nordlund 1997[ 1999), dense cores formed 
by shocks within IVlGs (Padoan et al.||2UUl) have an en¬ 
hanced magnetic-field strength on average. Furthermore, 
cores are topologically the ultimate zero-dimensional des¬ 
tination of a fluid element undergoing compression, as 
they can be viewed as the intersection of filaments that 
are formed by the intersection of postshock sheets. Much 
of the flow turbulent energy is dissipated by the time 
it ‘stagnates’ into a core. Due to this drop in turbu¬ 
lent energy, together with the increase in magnetic-field 
strength and density, the turbulence inside dense cores 


Fig. 25. — Alfvenic rm s M ach number versus cloud mass for the 
same clouds as in Figure !^ computed with the cloud mean mag¬ 
netic field (empty circles), or the cloud rms magnetic field (filled 
cirlces). 

is trans-Alfvenic and trans-sonic. In summary, super- 
Alfvenic and supersonic MG turbulence is the natural 
consequence of large-scale trans-Alfvenic trans-sonic tur¬ 
bulence and also the natural origin of small-scale trans- 
Alfvenic trans-sonic turbulence in prestellar cores. 

The small-scale enhancement of the magnetic field 
within MCs is partly illustrated in Figure by the val¬ 
ues of the rms field in the clouds that is approximately a 
factor of two larger than the mean field in the most mas¬ 
sive clouds. As a more direct demonstration of the super- 
Alfvenic nature of MG turbulence. Figure ^E\ shows the 
cloud rms Alfvenic Mach number versus the cloud mass. 
The Mach number is computed as the ratio of the cloud 
rms velocity and the cloud Alfven velocity, where the 
latter is computed either with the mean magnetic field 
(empty circles) or with the rms magnetic field (filled cir¬ 
cles), and using the mean density sampled by the tracer 
particles. Nearly all clouds with mass larger than 10^ 
Mq are super-Alfvenic, even considering their amplified 
field strength. For the 41 GMGs with masses larger than 
10^ Mq, the average Alfvenic Mach number is 8.3 with 
respect to the mean field, and 3.9 with respect to the rms 
field. 

The super-Alfvenic nature of the turbulence in the 
clouds from our simulation is consistent with the observa¬ 
tional evidence. Based on the comparison betwe en simu¬ 
lations of MHD turbulenc e and MG observations, |Padoan| 


& Nordlund (1997 1999) suggested that MG turbulence 
was better characterized by supersonic turbulent flows 
with ^ 1 than flows with A4a ~ 1- This result 


tions (Padoan et al. 

2004a|) and synthetic Zeeman split- 

ting measurements ( 

Lunttila et al.| 

2008 

2009). Taking 

advantage of the anisotropy of IVIHD tur 

bulence, |Heyer| 


demonstrated that the densest regions 
IG complex are characterized by super- 


& Brunt (|2012j) 
of the T'aurus 1\ 

Alfvenic turbulence, while in low density regions the mo¬ 
tions are sub or trans-Alfvenic, also consistent with the 
picture from our simulation, where MGs are formed by 
large-scale trans-Alfvenic turbulence, and thus fed pref¬ 
erentially by motions along magnetic field lines, as dis¬ 


cusse d above (jNordlund & Padoan||2003| jPadoan et al. 

mol) 
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Fig. 26. — B — n relation obtained from averaging the B — n 
relations of all 1563 clouds in our sample. The circles correspond 
to the mean value of uh and B in each density bin, while the error 
bars mark the values two standard deviations above and below the 
mean, to illustrate the scatter in the B — n relation. The B — n 
relation of an individual cloud is computed using the values of 
B and tt-h of all the tracer particles in that cloud. Because the 
clouds are selected as regions with density nn > 100 cm“^, no 
tracer particles belonging to a cloud has density lower than that. 
To extend the relation to densities tt-h < 100 cm“^, we define a 
cloud volume delimited by the smallest and largest coordinates of 
the tracer particles of that cloud, and include the values of B and 
riH of each cell of that volume to compute the relation. The two 
solid lines are least-square fits for densities nn < 10^ cm“^ and 
riH > 10^ cm“^, giving B ~ and B ~ respectively. 


To further characterize the cloud turbulence, we have 
computed the B — n relation inside all the clouds of our 
sample, using again the values of B and n of the tracer 
particles inside the clouds. We divide the density in 
logarithmic bins, and compute the mean magnetic field 
strength and its standard deviation in each density bin. 
We then average these values among all the clouds, us¬ 
ing weights proportional to the number of tracer particles 
in the density bins of each cloud. Figure shows this 
mean B — n relation. We also illustrate the large scatter 
by plotting error bars that correspond to twice the stan¬ 
dard deviation above and below the mean values. The 
two solid lines are least-square fits for nn < 10^ cm“^ 
B 


,0.13 


, and for nn > cm B 


’ n 


0.29 


The stronger dependence of the magnetic-field strength 
on density at nu > 10^ cm“^ than at lower density is 


qnalitativel y consistent with the observations (Crutcher 


et al. 


2010). T he slope we derive is m uch smaller than 


that derived by Crutcher et al. (2010) at high densities, 
B ^ However, their slope does not refer to the 

mean magnetic field at a given density, but to its max¬ 
imum value. Considering the large number of measured 
upper limits well below such maximum values, the de¬ 
pendence of the mean B on density could be significantly 
shallower than the estimated slope of the upper envelope 
of the B — n relation. Furth ermore, the Bayesian analy¬ 
sis by Crutcher et al. (2010) assumes a uniform distribu¬ 
tion of the magnetic lield strength, while this distribu¬ 


tion is exponentia l in super-Alfvenic turbulence (Padoan 
fc Nordlnnd|1999j ) (we have verified it is exponential also 


in our clouds). Finally, and most importantly, the B — n 
relation in Figure |26| is not computed for a selection of 
dense cores, as in the^observations, but using every single 
tracer particle in the cloud, so it should not be compared 


quantitatively to the observational B — n relation. Such 
a comparison would require synthetic Z eeman observa¬ 
tions of a se lection of dense cores, as in [Lunttila et al.| 
([20081 |2009|) . It would also require higher spatial resolu- 
tion, because most of the observed cores with the largest 
detected magnetic-field strengths, at densities of order 
10^ — 10^ cm“^, have sizes substantially smaller than the 
spatial resolution of our simulation. The higher resolu¬ 
tion would also allow to better resolve the dyna mo am¬ 
plification in dense cores ( Federrath et al.||2011b ), which 
would tend to increase the slope of the B — n relation. 

10. COMPARISON WITH OUTER-GALAXY MCS 

To further test our results, we carry out a comparison 
of the properties of our MCs with those of observed MCs. 
This is a preliminary approach based on the derivation of 
projected quantities, such as column density, equivalent 
radius and line-of-sight velocity dispersion. Follow-up 
studies with synthetic observations taking into consider¬ 
ation chemistry and radiative transfer are also needed. 

Our observational sample of choice for this comparison 
is the MC catalog by Heyer et al. (2001), extracted from 
a de composition of the FCRAO Outer Galaxy Sur¬ 
vey (Heyer et al.|1998). Besides the large dynamic range 
of the survey, its mam advantage is that for the Outer 
Galaxy there is no blending of emission from separate 
MGs along the line of sight, or at least the problem is 
strongly mitigated compared with the Inner Galaxy. As 
a result, a very large number of clouds can be reliably 
selected over a broad range of cloud masses and sizes. 
The catalog contains a total number of 10,156 objects, 
up to a mass of approximately 8 x 10 ^ Mq and a size of 
45 pc. It is estimated to be complete down to a mass of 
approximately 600 Mq and a cloud size of 3 pc. 

Inner-Galaxy MC catalogs are far less reliable and com¬ 
plete because of velocity blending, so they are not suit¬ 
able for the comparison we pursue her e. For example, 
the recent catalog of Inner-Gala xy MCs (Rathborne et al._ 
2009 Roman-Duval et al.|2009) extract^ from the UB- 

FGRAO Galactic Ring Survey ( Jackson et al.|2006 ) con¬ 

tains objects between 1 and lO^Tl07buFIi_^itmiated 
to be comp lete only above 4 x 10^ Mq (Roman-Duval 
et al. 2009). We suspect a more realistic completeness 
limit may be 2 x 10^ Mq, because t he mass distribution 


is a p ower law only above that mass (Roman-Duval et al. 
2009), which corresponds to a size limit of approximately 
20 pc, judging from the mass-size relation and from the 
lack of a power law in the size distribution below 20 pc. 
This severe incompleteness suggests that the cloud sur¬ 
face density may be overestimated by a large factor. The 
MC mass distribution is expected to be a power law down 
to small masses, so the abrupt departure from a power 
law below 2 x 10^ Mq indicates that much of the missing 
mass from smaller clouds is incorrectly assigned to larger 
ones due to velocity blending. The failure to select indi¬ 
vidual three-dimensional clouds is also demonstrated by 
the absence of a velocity-size correlation and, possibly, by 
the extremely low values of the cloud virial parameters 
(the distribution peaks at o^vir ~ 0.2), which would im¬ 
ply a larger star formation rate and a stronger signature 
of global collapse than observed. In summary, the dif¬ 
ferences between Galactic-Ring and Outer-Galaxy MCs 
may not be as large as often assumed. Thus, although 
our comparison is primarily with Outer-Galaxy MCs, the 
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Fig. 27. — Probability distribution of cloud mass from the sam¬ 
ple of 1,547 clouds of our highest-resolution catalog (512^ cells, or 
0.49 pc) and with the threshold density that best matches the ob¬ 
served mass-size relation, nH,min = 200 cm“^ (shaded histogram). 
The dashed line is the result of a least-square fit yielding a slope of 
—0.88 ± 0.06. The unshaded, solid-line histogram shows the mass 
distribution from the observational sample of 3,228 Outer-Galaxy 
MCs_(see text), selected from the larger sample in |Heyer et al.| 
(|2001|). The dotted-line histogram is the mass distribution ot the 
observational sample, excluding clouds with radius Re < 3.1 pc, 
the size-completeness limit of the Outer-Galaxy survey. It shows 
a contribution to the mass distribution by clouds below this com¬ 
pleteness limit up to a mass of approximately 1,000 Mq. The least- 
square-fit to the observational mass distribution above 1.5 X 10^ 
M 0 has a slope of —0.99 ±0.09 and is shown by the dashed-dotted 
line. 


results may be applicable to Galactic clouds in general. 
As in Heyer et al. (2001), we consider only the sub¬ 


set of 3,y01 clouds with circular velocities Vc < —20 km 
s“^, because of kinematic distance accuracy. We further 
select clouds with mass Md > 100 Mq, as that is the 
mass limit for our numerical cloud catalogs, resulting in 
a total sample of 3,228 Outer-Galaxy MGs. Given the 
distances to the clouds and the angular resolution of the 
survey, the spatial resolution varies between 0.4 pc and 
3.8 pc. Therefore, the cloud extraction of our highest 
resolution catalog with dx = 0.48 pc matches well the 
highest resolution in the observations. The main limita¬ 
tion of the survey is the velocity resolution, only slightly 
better than 1 km s“^, which, combined with the measure¬ 
ment of line width based on the equivalent width instead 
of the antenna-temperature-averaged velocity dispersion, 
results in a minimum velocity dispersion of clouds of ap¬ 
proximately 0.5 km s“^. However, we show that the data 
can be used to test both the slope and the normalization 
of the Larson velocity-size relation from the simulation 
despite this low velocity resolution. 

A explained in §3, we illustrate this comparison using 
clouds selected from 7 snapshots from the last 6 Myr 
of the simulation. However, we have verified that all 
the observational MG properties discussed in this section 
are essentially the same when derived from a catalog of 
clouds selected from the last 6 Myr prior to the inclusion 
of gravity. This confirms that global MG properties are 
primarily the result of SN-driven turbulence, with little 
modification due to self-gravity, apart from the slight in¬ 
crease in mean cloud density shown in §7, and an increase 
in the mass of the largest cloud. 

Of the 12 cloud catalogs described in §3, we choose 


Fig. 28. — Exponents of the power-law fits to the probability 
distributions of cloud masses (lower plots) and cloud sizes (up¬ 
per plots) for all twelve catalogs of MGs selected from the simula¬ 
tion. The exponents are plotted as a function of threshold density, 
^H,min 5 cind for three different values of spatial resolution of cloud 
selection. 


the one with the highest-resolution (512^ cells, or 0.49 
pc) and with the threshold density that best matches 
the observed mass-size relation, nH,min = 200 cm“^, for 
all the plots in this section, and we compute velocity dis¬ 
persions using the tracer particles, thus taking advantage 
of the highest resolution of the simulation, dx = 0.24 pc, 
due to the large number density of tracer particles within 
dense clouds. This catalog contains 1,547 objects. 


10.1. Mass Distribution 

Figure ^7\ shows the mass distribution of our clouds 
(shaded histogram). The histogram is well approximated 
by a power law, with a slope of —0.88 ± 0.06 in the ap¬ 
proximate mass range 200 — 10^ Mq (dashed line). The 
figure also shows the mass distribution of the MGs from 
the observational sample that is also nicely fit by a power 
law, with a slope of —0.91 ±0.09 in the approximat e mass 
range 1.5 x 1 0^ — 2 x 10^ Mq (dashed-dotted line). Heyer 


et al. 


(2001) derived a slope of —0.80 ± 0.03 by includ¬ 


ing all clouds down to the completeness limit of 600 Mq . 
Our slope is a bit steeper because we are a bit more 
conservative on the completeness limit. The dotted-line 
histogram in Figure shows the mass distribution for a 
sub-sample where we include only clouds above t he size 
comp l etenes s limit of 3.1 pc, the value derived by |Heyer| 
et al. (2001). The comparison with the histogram of the 
full sample shows that some MGs with sizes below the 
size completeness limit are found with masses up to ap¬ 
proximately 1,000 Mq, so we consider the catalog to be 
complete only above that mass value. 

The mass distribution of our clouds is consistent with 
that of the MGs from the Outer Galaxy Survey. Further¬ 
more, this is true for the mass distribution from all nu¬ 
merical cloud catalogs we have compiled, independent of 
the numerical resolution and threshold density of cloud 
selection. As shown in Figure the exponent of the 
power-law fit to the mass distribution has a weak depen¬ 
dence on resolution and threshold density and, within 
the 1-cr error bars shown in the plot, it is consistent with 
the observational exponent in all cases. 

The largest cloud in our nH,min = 200 cm“^ catalog 
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Fig. 29. — Probability distribution of cloud size for the same 
cloud catalog from the simulation (shaded histogram) and the same 
obs ervational sample (unshaded, solid-line histogram) as in Figure 
|27| The dotted-line histogram is the size distribution of the obser¬ 
vational sample, excluding clouds with mass M^i < 600 Mq, the 
mass-completeness limit of the Outer-Galaxy survey. The dashed 
line is a fit to the tail of the histogram from the simulation, with 
slope —2.5 ± 0.2, and the dashed-dotted line a fit to the observa¬ 
tional size distribution, with slope —2.3 ± 0.3. 


has a mass of 1.3 x 10^ Mq, a few times smaller than 
the largest MC in the observational sample. However, 
Figure shows that this maximum mass is consistent 
with the largest mass expected from our sample size and 
the slope of the mass distribution. If we simulated a 
region larger than 250 pc, and thus collected a much 
larger cloud sample, we would likely derive a power-law 
mass distribution extended to larger masses. The obser¬ 
vations are consistent with a power-law mass distribution 
for clouds more massive than those in the Outer Galaxy 
Surve y. For example, in the i r analysis of the MC sam- 
ple by Solomon et al. (1987), Williams & McKee (1997) 
found a comparable slope of —0.81 ± 0.14 in the mass 
range 3 x 10^ — 5.6 x 10^ Mq. Although more uncertain, 
th e slope they obtained from the analysis of the sample 
by |Scoville et al.| (|1987|), —0.67 ± 0.25, is also consistent 
with the slope of tne Outer-Galaxy MCs a t lower masses. 
More recently, Roman-Duval et al. (2010) found a slope 
of —0.64 ± 0.25 in the mass range 4 x 10 -10^ Mq from 
the analysis of the Galactic Ring Survey (a more conser¬ 
vative completeness limit of 8 x 10^ Mq gives a slope of 
-0.86 ±0.25). 


10.2. Size Distribution 

Following |Heyer et al.| (|2001|) and most observational 
works, as a measure of a cloud's size we adopt the equiv¬ 
alent radius. Re = (Aci/tt)^’^, where is the cloud 
projected area. The probability distribution of Re for 
the clouds from the simulation is shown in Figure 
(shaded histogram). It is well approximated by a power 
law with a slope of —2.5 ± 0.2 in the approximate range 
of 2 — 15 pc. Within the uncertainty, this is consistent 
with the slope of —2.3 ± 0.3 of the observational sam¬ 
ple in the approximate equivalent-radius range of 5 — 50 
pc. Furthermore, Figure 28\ shows that the size distribu¬ 
tions of clouds selected from the simulation with differ¬ 
ent threshold density and resolution are also consistent 
with the observations, within the 1-cr uncertainty (ex¬ 
cept for the catalog with nH,min = 200 cm“^ and the 


Fig. 30.— Equivalent radius versus the three-dimensional cloud 
rad ius for a ll the clouds in the same simulation sample as in Figures 

Ell and EH 


lowest-resolution). As in the case of the mass distri¬ 
bution, we have been slightly more conservative in the 
estimation of the size completeness limit, based on evi¬ 
dence that around the value of 3.1 pc , the completeness 
limit estimated by Heyer et al. (2001), we still find some 
contribution from clouds with mass below the mass com¬ 
pleteness limit of 600 Mq, as illustrated by the dotted 
histogram i n Figure M As a result, we find the same 
slope as in |Heyer et al.| (|2001|) , but with a three times 
larger uncertainty. T'he same powe r law seems to apply 
to even larger clouds. For example, Sanders et al. (1985) 
find a slope of —2.3 ±0.2 for a sample of 80 clouds in the 
approximate size range of 20 — 80 pc. 

Because we have previously computed a three dimen¬ 
sional cloud radius, Rcb from the simulation, we can test 
the relation between the observable radius, Re^ and the 
three-dimensional one. The comparison is shown in Fig- 
ure[^ where the dashed line corresponds to Re = Rc\- 
Re is smaller than Re\ for most clouds with Rd ^ 2 
pc, and larger than Rc\ for most clouds with Rc\ < 2 
pc. The average ratio is Re/Rc\ = 0.87 and increases 
towards smaller radii. As a consequence, the probability 
distribution of Rc\ is slightly shallower than that of Re. 


10.3. Velocity-Size Relation 

The velocity-size relation of MGs determines the nor¬ 
malization of the velocity scaling inside individual clouds. 
In §6.2, we have shown that the energy from SN ex¬ 
plosions sets a turbulent cascade inside individual MGs 
that follows the usual velocity scaling of supersonic tur¬ 
bulence, but we have not discussed the normalization 
of the velocity structure functions of individual clouds. 
To address the velocity normalization, we compute the 
internal rms velocity of our clouds based on the veloc¬ 
ity of their tracer particles. Being derived from tracer 
particles, this rms velocity is mass-weighted, which is 
a reasonable approximation when comparing it with the 
antenna-temperature-weighted rms velocity from MG ob¬ 
servations. As in the observations, we compute the one¬ 
dimensional (line-of-sight) rms velocity, in the direction 
perpendicular to the plane (plane of the sky) where we 
measure the equivalent radius. 

Figure [3T] shows the velocity-size relation of the clouds 
from the simulation (empty circles) and from the obser- 
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Fig. 31. — One-dimensional rms velocity versus size of clouds 
selected from the simulation (empty circles) and from the Outer 
Galaxy Survey ( fille d ci rcles). The cloud samples are the same 
as in Figures [27| and but the observations are shown only for 
i^e > 4 pc, because for smaller cloud sizes the lower envelope of 
the velocity-size relation is not resolved by the observation (the 
minimum value of ctv that can be detected is ~ 0.5 km s“^). The 
thin solid line shows the mean values of ctv in logarithmic bins of 
Re: and the thick solid line is a fit to those values, giving a slope of 
0.39 ± 0.07. The dashed line is the fit to the binned data from the 
observations, giving a slope of 0.48 ± 0.06, and exactly the same 
normalization as the simulation at Re ~ 10 pc. 


vational sample (filled circles). We have excluded all the 
observed clouds with Rq < 4 pc, because the observations 
cannot detect velocity dispersions smaller than approx¬ 
imately 0.5 km s“^, due to the low velocity resolution. 
The lower envelope of the velocity-size relation of both 
the simulation and the observations decreases with de¬ 
creasing cloud size, and reaches the value of 0.5 km s“^ 
at approximately 4 pc. Thus, the velocity dispersion of 
a fraction of the clouds smaller than 4 pc may be signifi¬ 
cantly overestimated, with that fraction growing towards 
smaller cloud radii, causing an artificial flattening of the 
velocity-size relation. As shown in Figure 6 of |Heyer| 


et al. (2001), the velocity-size relation for the full sample 
IS essentially flat below 4 pc. 

Values of from the simulation are not to be trusted 
for the smallest cloud sizes, Rq < 2 pc, because of the in¬ 
creasing effect of numerical dissipation towards smaller 
scales. The thin solid line in Figure shows the av¬ 
erage values of in logarithmic bins of Rq. While it 
is nicely fit by a power law for Rq > 2 pc, it clearly 
drops at smaller cloud sizes. This is consistent with the 
cloud structure functions shown in Figure where the 
numerical dissipation starts to become important below 
approximately 2 pc as well. 

Despite these limitations imposed by the low velocity 
resolution of the observations and the numerical dissi¬ 
pation in the simulation, we still have a sufficient range 
in Rq where the simulation and the observations can be 
compared. Both the upper and the lower envelopes of the 
velocity-size relation are very similar in the two cases. 
Furthermore, the velocity normalization is nearly identi¬ 
cal. The thick solid and dashed lines in Figure [sT] show 
the least square fits of the average values of in loga¬ 
rithmic bins of Rq of the simulation (for Rq > 2 pc) and 
of the observations (for Rq > 4 pc), respectively. From 



10 100 

[Mo pcq 


Fig. 32.— Normalization of the velocity-size relation versus col¬ 
umn density. We normalize the velocity dispersi on with R^'^, in¬ 
stead of the shallowe r slopes derived in Figure to reproduce 
the plot in Figu re 7 of|Heyer et al.|(|2009|). Besides the data points 
from Figure!^ (this time including Outer-Galaxy MGs with i7e < 4 
pck we also show a sample of MGs from the Galactic Ring Survey 

i Roman-Duval et al.||2010|). The column density of the Galactic- 
tmg IViUs IS on the average 10 times larger than that of the Outer- 
Galaxy MGs, yet the velocity normalization is essentially the same. 

the simulation we get 

CTv = (1.34 ± 0.04) kIns“\i^e/10pc)°•^^=^°■°^ (19) 

and from the observations: 


(1.34 ± 0.06) kms“^(i?e/10pc) 


\0.48±0.06 


( 20 ) 


SO the velocity normalization at Rq = 10 pc is indistin¬ 
guishable in the two cases. This agreement between the 
simulation and the observations in the slope, total scatter 
and normalization of the velocity-size relation is strong 
evidence that SN driving alone can be responsible for the 
turbulence observed in MCs. 

The universality of the MC ve l ocity normalization has 
been questioned by Heyer et al. (2009), claiming that it 
depends on column density, and thus that the velocity- 
size relation is controlled by gravity rather than being 
a natural consequence of the ISM turbulence. To fur¬ 
ther confirm the agreement between the simulation and 
the observations, we show the velocity normalization 
as a function of column density in Figure 32 (we plot 
O’v as in|Heyer et al.| (|2009|), even if the slope of the 
velocity-size relation is actually smaller than 0.5). Be¬ 
cause the range of cloud column densities is similar in 
the two cases, we should not expect a different normal¬ 
ization even if it depended on surface density. Figure [32] 
shows a good overlap between our clouds and the obser¬ 
vations. We also plot the values for the Galactic-Ring 
clouds in |Roman-Duval et al.| (|2010|) that have an av¬ 
erage surface density an order of magnitude larger than 
the Outer-Galaxy clouds (notice that the difference in 
surface density is only a factor of five for equal cloud 
mass or size, and could have been overestimated by a fac¬ 
tor of two or three, as explained in the opening of §10). 
The figure shows that there is no difference in the ve¬ 
locity normalization of Galactic-Ring and Outer-Galaxy 
MGs, despite the difference in surface density. Thus, we 
conclude that the normalization of the velocity-size rela¬ 
tion of the MGs in our sample is consistent with being 
controlled by SN-driven turbulence, rather than by the 
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Fig. 33. — Cloud rms velocity versus cloud size for the same 
numerical cloud sample as in the previous figures, but with both 
rms velocity and cloud size measured in three-dimensional space, 
CTv = ctv^sd/'s/S and The thin solid line shows the average 
in logarithmic intervals of Rd, and the dashed line is the fit to the 
binned values, with a slope 0.37 ±0.02 for cloud sizes Rd > 1.6 pc. 


clouds self-gravity. This result is in contradiction to the 
claim that the v elocity nornializat ion of MCs scales with 
surface density dHeyer et al.|2QQ9|), ba sed on clouds from 
the sample by polomon et al.| (|i9b7p, analyzed within 
rectangular maps of ditterent sizes, rather than a fixed 
antenna-temperature threshold. 

Besides being useful to derive the velocity normaliza¬ 
tion, the observed velocity-size relation may also pro¬ 
vide a rather accurate estimate of the slope of the sec¬ 
ond order velocity structure function. In Figure we 
show the relation between the velocity dispersion derived 


from the three-dimensional one, cTv^sd/ and the three- 
dimensional cloud radius, Rc\. The thick solid line is the 
least square fit to the average values of ay in logarithmic 
bins of Rc\ shown by the thin solid line. The fits to the 
three-dimensional velocity-size relation for Rc\ > 1.6 pc 
(where the relation is well approximated by a power law) 
gives 


CTv = (1.19 ± 0.04) kms“^(i?ci/10pc)°-3^=^° °2. (21) 


This relation is consistent with its two-dimensional coun¬ 
terpart (apart from the lower normalization due to the 
fact that Rc\ > Re on average), so the observable rela¬ 
tion can be considered as a good estimate of the intrinsic 
three-dimensional one. Furthermore, the slope of the re¬ 
lation (21) is also consistent with the slope of the second 
order structure function, C(3)/2 ~ 0.37, of the clouds 
from the simulation. Thus, we conclude that the ob¬ 
served velocity-size relation provides an estimate of the 
velocity scaling of MC turbulence, as long as it is based 
on MCs with reliable distance measurements and suffi¬ 
cient velocity resolution to detect the lower envelope of 
the relation. 


The velocity-size relation (21) implies the following ex¬ 


pression for the dynamical time of MCs as a function of 
the cloud three-dimensional radius, using the definition 
(13) of dynamical time adopted in §8: 


idyn = 4.8Myr (i?ci/10pc) 


0.63 


( 22 ) 


Using the result of §8 that {tine/tdyn) ~ 4, our velocity- 
size relation implies that our largest MCs with sizes in 


Fig. 34. — Mass versus size for t he s ame numerical and obser¬ 
vational cloud samples as in Figure |31| but including all observed 
MCs with mass larger than 100 Mq. The observational points 
have been shifted to the left by dividing the observed values of Re 
by 10, to avoid the overlap with the numerical data points. The 
thick solid and dashed lines are fit to the binned data of the simu¬ 
lation and the observations respectively, for Re > 2 pc. The thin, 
long-dashed lines show two values of constant surface density. 

the range Rc\ ^ 10 — 30 pc have lifetimes in the range 
tiife ^19 — 38 Myr. 


10.4. Mass-Size Relation 

The mass-size relation is plotted in Figure [34l this time 
including all observed clouds above 100 Mq. The values 
of Re of the observed MCs have been divided by 10, to 
avoid overlap with the clouds from the simulation. Be¬ 
cause of the imposed limit on the minimum cloud mass, 
the data is binned and fit only for Re > 2 pc for both the 
simulation and the observations. The resulting fits are 

Mci = (9.6 ± 0.3) X lO^M© (i?e/10 (23) 

for the simulation, and 

Mci = (10.9 ± 0.7) X lO^M© (i?e/10pc)2'49±o-07^ (24) 


for the observations, so the slope of the mass-size relation 
from the simulation is consistent with the observations. 
The normalization of the mass-size relation depends on 
the threshold antenna temperature of the observational 
sample and the threshold density of the numerical sam¬ 
ple. Of our MC catalogs described in §3, the ones with 
^H,min = 200 cm“^ havc the mass-size normali zation 
closest to tha t of the Outer-Galaxy MC sample by |Heye^ 


et al. (2001). This is the reason why all plots in this 
section are based on the highest-resolution catalog with 
that value of threshold density. 

The total scatter in the relation is also similar between 
the observations and the simulation, if we account for the 
facts that the observational sample size is approximately 
twice as large as the numerical one, and that we have 
not added any observational uncertainty to the masses 
and sizes of the clouds from our simulation. Because 
of the similarity in both the slope and the total scat¬ 
ter, we can conclude that the mass-size relation resulting 
from SN-driven turbulence is consistent with that of real 
MCs from the Outer Galaxy Survey. A similar mass-size 
relation (though with a five times larger norma lization) 


was also derived f rom the Galactic Ring Survey (Roman- 
Dnval et al.|2QlQ ) and is implied by previous estimates of 
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Fig. 35. — V irial parameter versus mass for the same clouds as 
in Figures [M| The dashed line shows a predicted upper envelope 
assuming that the largest velocity dispersion is ~ 2.5 km s“^ (inde¬ 
pendent of cloud mas s) and adopting the fit to the mass-size rela¬ 
tion from Figures [Ml The dashed-dotted line shows the predicted 
lower envelope, also based on the mass-size relation and assuming 
a minimum velocity of ~ 0.5 km s“^, as in the observations. 



cloud fractal dimensions from various observational sur- 


veys (e.g;. Elmegreen & Falgarone 1996 Sanchez et al. 
2QQ7| ) and from simulations of randomly driven tur bu¬ 


lence (e.g. [Kritsuk et al.||2007t jEederratli et al.||2009). 

Combining the mass-size relatio n p3| ) with the dynam- 
ical time expression in equation ([22]h and adopting the 
average value derived in §10.2 for the ratio between the 
two definitions of cloud radius, Rq = 0.87i?cb we obtain 
the following expression for the average dynamical time 
of MCs as a function of cloud mass, 


idyn = 5.OMyr(Mci/lOm0f25, (25) 

hence the expression ( |18| ) for the average cloud lifetime 
as a function of cloud mass anticipated in §8. 


10.5. Virial Parameter 

In §7, the virial parameter is computed with the three- 
dimensio nal radius, Rc ] . Her e, as in most MC studies 
including Heyer et al. (2009), we define an observable 
version of the virial parameter using the equivalent ra¬ 
dius, i?e, and thus refer to it as (avir,e- The dependence of 
(avir,e on Tfci is fully determined by the velocity-size and 
mass-size relations presented above. Nevertheless, it is 
presented here as an alternative view of the comparison 
of the simulation with the observations, and to suggest 
an empirical calibration of the virial parameter of MCs. 
The values of avir,e are plotted versus in Figure 


35 


va^ 


for both the simulation (empty circles) and the obser- 
lons (filled circles). The most striking feature of this 
plot is the very large scatter in the values of (avir,e at a 
fixed cloud mass, growing with decreasing cloud mass, 
which is a direct consequence of the large scatter in the 
velocity-size relation. As in that relation, the lower en¬ 
velope for the observational data is limited by the small¬ 
est value of cTv to which the observations are sensitive, 
(Tv ~ 0.5 km s“^. The dashed-dotted line in Figure 35 
shows the minimum value of (avir,e as a function of Md, 
computed by setting Rq as a function of Md using the 
average mass-size relation from the previous section and 
setting (Tv = 0.5 km s“^. 


The virial parameter would be independent of mass if 

the velocity size relation were Re and the mass- 

size relation were Md ^ i?g, the often assumed form 
of the Larson relations. However, the mass-size rela¬ 
tion is such that the cloud surface density grows with 
mass, as shown in the previous section, causing the de¬ 
crease of avir,e with increasing Md seen in Figure]^ only 
partly mitigated by the exponent of the velocity-size re¬ 
lation being a bit smaller than 0.5. The upper envelope 
in Figure is even steeper than the average decrease 
of virial parameter with mass, as a consequence of the 
nearly flat upper envelope of the velocity-size relation. 
The dashed line in Figure |35| shows the expected upper 
envelope based on the mass-size relation and a maximum 
rms velocity of 2.5 km s“^ that is representative of some 
of the largest values in both the simulation and the ob¬ 
servations. 

Once we account for the minimum velocity disper¬ 
sion in the Outer-Galaxy Survey, the relation of (avir,e 
with Md and its scatter for the clouds from the simu¬ 
lation are consistent with those for the observed Outer- 
Galaxy MCs, as expected from the agreement found in 
the velocity-size and mass-size relations. This agreement 
suggests the possibility of an empirical calibration of the 
observed values of the virial parameter based on the 
results of our simulation. We have shown in §7 that 
the virial parameter computed with the radius Rc\ is 
Q^vir ~ 1.2 (2Fk/Fg). We have also shown in this section 
that, on the average, Rd ~ i?e/0.87, thus we derive this 
useful result for the relation between the observed virial 
parameter based on the equivalent cloud radius and the 
energy ratio: 

2Ev 

« 0.96Q!vir,e (26) 


(assuming negligible saturation of the observed lines, 
such that the rms velocity c an be assumed to be ap- 
proximately mass weighted). Bertoldi & McKee (1992) 
modeled extensively the coefficient of the virial param- 
eter of clumps, depending on the mass profile and the 
ellipticity of the clumps. Given their complex structure, 
MCs are not described by radial density profiles and el- 
lipticity parameters as easily as compact clumps, so an 
empirical calibration based on simulations as proposed 
here is needed. 


11. OVERVIEW OF RESULTS AND CONCLUSIONS 

The main goal of this work is to test if ISM turbulence 
driven only by SN explosions can explain the turbulence 
observed within MCs. We have addressed this question 
with an AMR simulation representing an ISM volume 
of (250 pc)^ and reaching a maximum resolution of 0.24 
pc, with refinement based on density, density gradients 
and pressure gradients to resolve individual SN remnants 
and their interaction with the dense gas. We have stud¬ 
ied the SN-driven turbulence over the whole volume and 
within individual clouds. We have compiled 12 differ¬ 
ent catalogs of MCs selected from the simulation using 
four different values of threshold density and three dif¬ 
ferent spatial resolution. The properties of these clouds 
have been studied using tracer particles, hence taking 
advantage of the full resolution of the simulation. First 
we have presented cloud properties based on particle po¬ 
sition, velocity and magnetic field values measured in 
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three-dimensional space; then we have carried out a com¬ 
parison with real MCs from the Outer Galaxy Survey by 
measuring projected quantities, such as the line-of-sight 
velocity dispersion and the equivalent radius. Our results 
are summarized in the following: 

1. Near equipartition of total energies in the whole 
volume results in a distinct energy separation in 
the dense gas. While the overall ISM turbulence 
is trans-Alfvenic and mildly supersonic, the tur¬ 
bulence in the dense gas is highly supersonic and 
super-Alfvenic. 

2. Approximately 11% of the total kinetic energy is 
transferred to the dense gas, even if most SN ex¬ 
plosions occur at low densities. The dense gas has 
an average velocity dispersion of 8.5 km s“^. 

3. During the rapid expansion of a SN remnant, the 

velocity power spectrum may briefly show strong 
features at different scales. During most of the time 
the power spectrum is statistically relaxed and de¬ 
velops a power-law inertial range that scales with 
wavenumber as E{k) oc between approxi¬ 

mately 20 and 60 pc. 

4. Unlike previous studies of compressible turbulence, 

the power spectrum of the compressive component 
of the velocity, Ec{k) oc is much steeper 

than that of the solenoidal component, Es{k) oc 
^-1.31 rppg baroclinic effect is the best candidate 
to explain this result, as previous simulations of 
supersonic turbulence adopted equations of state 
where such effect was absent. 

5. SN driving is not purely compressive. The curl of 
the forcing is the baroclinic term that is compa¬ 
rable to or larger than the divergence of the forc¬ 
ing (except in the unrealistic cases of a uniform- 
density medium or a barotropic equation of state). 
As a result, the power in solenoidal modes ex¬ 
ceeds that in compressive modes almost at any 
time and any wavenumbers (at the scale of 20 pc, 
Ec/Es ^ 0.2). Thus, isothermal simulations of tur¬ 
bulent fragmentation based on random solenoidal 
driving are a much better approximation of MC 
turbulence than isothermal simulations with purely 
compressive driving. 

6. The time-averaged energy-injection scale of SN- 
driven turbulence is approximately 70 pc with our 
SN rate (may be larger with a smaller SN rate, or 
vice-versa), and oscillates in time between 50 and 
100 pc. 

7. The scaling exponents of the first and second order 
structure functions of velocity, relative to the third 
order one, in SN-driven turbulence are consistent 
with those found in supersonic turbulence driven by 
an idealized, random volume force, which supports 
the validity of turbulent fragmentation and star for¬ 
mation studies where the ISM turbulent cascade 
within MCs was modeled with such an idealized 
driving. 


8. Based on a new scheme to compute density- 
weighted velocity structure functions, we obtain a 
third order exponent close to unity, as in the in¬ 
compressible Kolmogorov’s ‘4/5 law’. Thanks to 
the AMR method and to this density-weighting 
scheme, the structure functions probe the inertial 
range of MC turbulence down to a scale of 2-3 pc, 
while previous studies of SN-driven turbulence did 
not resolve an inertial range and only addressed the 
relative scaling. 

9. The scaling of the velocity structure functions 
within individual MCs selected from the simula¬ 
tion is generally consistent with the scaling derived 
from MC observations. However, deviations are 
found for MCs directly affected by recent SN rem¬ 
nants. 

10. The ratio of cloud virial parameter and kinetic to 
gravitational energy ratio is (avir/(2^k/^g) = 1-2, 
independent of energy ratio and mass (see point 15 
for the observable virial parameter). This struc¬ 
tural property of MCs is not significantly affected 
by self-gravity during the duration of our simu¬ 
lation, where self-gravity is included in the final 
11 Myr, corresponding to about two average cloud 
free-fall times. The ratio becomes very large only 
in a small fraction of clouds with o^vir ^10, where 
self-gravity causes the collapse of dense cores. Even 
in these star-forming clouds, there is no evidence 
of global cloud collapse. Self-gravity only causes a 
slight shift towards larger densities of the probabil¬ 
ity distribution of cloud densities. 

11. The formation and dispersion times of MCs are of 
the order of two cloud dynamical times. Equiva¬ 
lently, the cloud lifetime, defined as the sum of for¬ 
mation and dispersion times, is approximately four 
cloud dynamical times. This is evidence indicating 
that SN-driven turbulence is responsible for cloud 
formation and dispersion, with little influence from 
self-gravity visible during the duration of our run. 
Euture work, with longer runs, is needed to deter¬ 
mine to what extent this remains true for longer 
periods of time. 

12. The clouds have a mean magnetic field enhanced 
only by a factor of two relative to the mean mag¬ 
netic field in the simulated volume The turbulence 
is super-Alfvenic for all clouds more massive than 
approximately 10^ Mq. 

13. The comparison with the MC sample from the 
Outer Galaxy Survey shows that clouds selected 
from the simulation have properties consistent with 
the observations, such as the mass and size distri¬ 
butions, the velocity-size and mass-size relations 
and the dependence of virial parameter on cloud 
mass. In our run, these properties, including the 
normalization of the velocity-size relation, are es¬ 
sentially the same for clouds selected either before 
or after the inclusion of gravity in the simulation; 
they are primarily the result of SN-driven turbu¬ 
lence, with only a minor contribution from self¬ 
gravity. 
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14. The normalization of the velocity-size relation does 
not depend on surface density in the simulation, 
nor in the observations. It is the same for MCs 
from the Outer Galaxy Survey as for those in the 
Galactic Ring Survey whose surface density is sig¬ 
nificantly higher. 

15. The simulation provides a calibration of the observ¬ 
able virial parameter, based on the equiva¬ 

lent radius, which allows a derivation of the en¬ 
ergy ratio from observational quantities, ^ 

0.96 (avir,e- 

Based on these results from the simulation and given 
its successful comparison with the Outer Galaxy Survey, 
we conclude that the SN-driven turbulence in our sim¬ 
ulation is consistent with the observed MG turbulence 
during the duration of our experiment. Although other 
energy sources are present in galaxies, and local radia¬ 
tive and mechanical feedbacks also play a role in the dis¬ 
persion of star-forming MGs, the origin, evolution and 
internal dynamics of MGs in our run are primarily the 
consequence of SN-driving, which is able to sustain tur¬ 
bulence at observed levels without help from those extra 
energy sources. 


We are grateful to Ghristoph Federrath and the anony¬ 
mous referee for providing useful comments and correc¬ 
tions. Gomputing resources for this work were provided 
by the NASA High-End Gomputing (HEG) Program 
through the NASA Advanced Supercomputing (NAS) 
Division at Ames Research Genter, by PRAGE through 
a Tier-0 award providing us access to the computing 
resource SuperMUG based in Germany at the Leibniz 
Supercomputing Genter, and by the Port d’Informacio 
Gientffica (PIG), Spain, maintained by a collaboration 
of the Institut de Efsica d’Altes Energies (lEAE) and the 
Gentro de Investigaciones Energeticas, Medioambientales 
y Tecnologicas (GIEMAT). PP acknowledges support by 
the ERG EP7-PEOPLE- 2010- RG grant PIRG07-GA- 
2010-261359 and by the Spanish MINEGO under project 
AYA2014-57134-P. TH is supported by a Sapere Aude 
Starting Grant from The Danish Gouncil for Independent 
Research. Research at Gentre for Star and Planet Eorma- 
tion was funded by the Danish National Research Eoun- 
dation and the University of Gopenhagen’s programme 
of excellence. 


REFERENCES 


Balsara, D. S., Kim, J., Mac Low, M.-M., &: Mathews, G. J. 2004, 
ApJ, 617, 339 

Benzi, R., Ciliberto, S., Tripiccione, R., Baudet, C., Massaioli, F., 
& Succi, S. 1993, Phys. Rev. E, 48, 29 
Bertoldi, F., & McKee, C. F. 1992, ApJ, 395, 140 
Boldyrev, S. 2002, ApJ, 569, 841 

Boldyrev, S., Nordlund, A., & Padoan, P. 2002, ApJ, 573, 678 
Bournaud, F., Elmegreen, B. G., Teyssier, R., Block, D. L., & 
Puerari, 1. 2010, MNRAS, 409, 1088 
Chevalier, R. A., &: Klein, R. I. 1978, ApJ, 219, 994 
Couch, S. M., Wheeler, J. C., & Milosavljevic, M. 2009, ApJ, 696, 
953 

Crutcher, R. M., Wandelt, B., Heiles, C., Falgarone, E., & 
Troland, T. H. 2010, ApJ, 725, 466 
de Avillez, M. A., & Breitschwerdt, D. 2005, A&A, 436, 585 
—. 2007, ApJ, 665, L35 
Dobbs, C. L. 2015, MNRAS, 447, 3390 
Dobbs, C. L., & Pringle, J. E. 2013, MNRAS, 432, 653 
Dubrulle, B. 1994, Phys. Rev. Lett., 73, 959 

Elmegreen, B. G., Elmegreen, D. M., &: Leitner, S. N. 2003, ApJ, 
590, 271 

Elmegreen, B. G., & Falgarone, E. 1996, ApJ, 471, 816 
Falkovich, G., Fouxon, L, & Oz, Y. 2010, Journal of Fluid 
Mechanics, 644, 465 

Faucher-Giguere, C.-A., Quataert, E., & Hopkins, P. F. 2013, 
MNRAS, 433, 1970 

Federrath, C. 2013, MNRAS, 436, 1245 

Federrath, C., Chabrier, G., Schober, J., Banerjee, R., Klessen, 

R. S., & Schleicher, D. R. G. 2011a, Physical Review Letters, 
107, 114504 

Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156 
—. 2013, ApJ, 763, 51 

Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79 
—. 2009, ApJ, 692, 364 

Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., & 
Klessen, R. S. 2011b, ApJ, 731, 62 
Ferland, G. J., Korista, K. T., Verner, D. A., et al. 1998, PASP, 
110, 761 

Fleck, Jr., R. C. 1996, ApJ, 458, 739 
Franco, J., & Cox, D. P. 1986, PASP, 98, 1076 
Gatto, A., Walch, S., Low, M.-M. M., Naab, T., Girichidis, P., 
Glover, S. C. O., Wiinsch, R., Klessen, R. S., Clark, P. C., 
Baczynski, C., Peters, T., Ostriker, J. P., Ibanez-Mejia, J. C., 

& Haid, S. 2015, MNRAS, 449, 1057 


Glover, S. C. O., & Clark, P. C. 2012, MNRAS, 421, 116 
Gnedin, N. Y., & Hollon, N. 2012, ApJS, 202, 13 
Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421 
Heiles, C., & Troland, T. H. 2005, ApJ, 624, 773 
Hennebelle, P., & Chabrier, G. 2008, ApJ, 684, 395 
—. 2011, ApJ, 743, L29 

Herant, M., & Woosley, S. E. 1994, ApJ, 425, 814 
Heyer, M., Krawczyk, C., Duval, J., &: Jackson, J. M. 2009, ApJ, 
699, 1092 

Heyer, M. H., Brunt, C., Snell, R. L., Howe, J. E., Schloerb, F. P., 
& Carpenter, J. M. 1998, ApJS, 115, 241 
Heyer, M. H., & Brunt, C. M. 2004, ApJ, 615, L45 
—. 2012, MNRAS, 420, 1562 

Heyer, M. H., Carpenter, J. M., & Snell, R. L. 2001, ApJ, 551, 852 
Heyer, M. H., & Terebey, S. 1998, ApJ, 502, 265 
Hopkins, P. F. 2012, MNRAS, 423, 2037 
lanjamasimanana, R., de Blok, W. J. G., Walter, F., Heald, 

G. H., Caldu-Primo, A., & Jarrett, T. H. 2015, AJ, 150, 47 
Iffrig, O., & Hennebelle, P. 2015, A&A, 576, A95 
Jackson, J. M., Rathborne, J. M., Shah, R. Y., Simon, R., Bania, 
T. M., Clemens, D. P., Chambers, E. T., Johnson, A. M., 
Dormody, M., Lavoie, R., & Heyer, M. H. 2006, ApJS, 163, 145 
Joung, M. K. R., &: Mac Low, M.-M. 2006, ApJ, 653, 1266 
Joung, M. R., Mac Low, M.-M., & Bryan, G. L. 2009, ApJ, 704, 
137 

Kifonidis, K., Plewa, T., Scheck, L., Janka, H.-T., &: Muller, E. 
2006, A&A, 453, 661 

Kim, C.-G., & Ostriker, E. C. 2015, ApJ, 802, 99 
Kolmogorov, A. N. 1941, Dokl. Akad. Nauk. SSSR, 30, 301 
Kowal, G., &: Lazarian, A. 2007, ApJ, 666, L69 
Kritsuk, A. G., Lee, C. T., & Norman, M. L. 2013, MNRAS, 436, 
3247 

Kritsuk, A. G., Norman, M. L., Padoan, P., &: Wagner, R. 2007, 
ApJ, 665, 416 

Kritsuk, A. G., Ustyugov, S. D., Norman, M. L., & Padoan, P. 
2010, in Astronomical Society of the Pacific Conference Series, 
Vol. 429, Numerical Modeling of Space Plasma Flows, 
Astronum-2009, ed. N. V. Pogorelov, E. Audit, & G. P. Zank, 
15 

Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250 
Larson, R. B. 1981, MNRAS, 194, 809 
Lazarian, A., & Pogosyan, D. 2000, ApJ, 537, 720 
Lehnert, M. D., Le Than, L., Nesvadba, N. P. H., van Driel, W., 
Boulanger, F., & Di Matteo, P. 2013, A&A, 555, A72 



Supernova Driving. 1. The Origin of Molecular Cloud Turbulence 


29 


Lunttila, T., Padoan, P., Juvela, M., & Nordlund, A. 2008, ApJ, 
686, L91 

—. 2009, ApJ, 702, L37 

Martizzi, D., Faucher-Giguere, C.-A., & Quataert, E. 2015, 
MNRAS, 450, 504 

Molina, F. Z., Glover, S. G. O., Federrath, G., & Klessen, R. S. 
2012, MNRAS, 423, 2680 

Monin, A. S., & laglom, A. M. 1975, Statistical fluid mechanics: 
Mechanics of turbulence. Volume 2 /revised and enlarged 
edition/ 

Neufeld, D. A., Lepp, S., & Melnick, G. J. 1995, ApJS, 100, 132 
Nordlund, A., & Padoan, P. 2003, in Lecture Notes in Physics, 
Berlin Springer Verlag, Vol. 614, Turbulence and Magnetic 
Fields in Astrophysics, ed. E. Falgarone & T. Passot, 271-298 
Ostriker, E. G., McKee, G. F., & Leroy, A. K. 2010, ApJ, 721, 975 
Ostriker, E. G., & Shetty, R. 2011, ApJ, 731, 41 
Padoan, R, Haugbplle, T., & Nordlund, A. 2012, ApJ, 759, L27 
—. 2014, ApJ, 797, 32 

Padoan, P., Jimenez, R., Juvela, M., & Nordlund, A. 2004a, ApJ, 
604, L49 

Padoan, P., Jimenez, R., Nordlund, A., & Boldyrev, S. 2004b, 
Physical Review Letters, 92, 191102 
Padoan, P., Juvela, M., Goodman, A. A., & Nordlund, A. 2001, 
ApJ, 553, 227 

Padoan, P., Juvela, M., Kritsuk, A., & Norman, M. L. 2006, ApJ, 
653, L125 

Padoan, P., Kritsuk, A. G., Lunttila, T., Juvela, M., Nordlund, 
A., Norman, M. L., & Ustyugov, S. D. 2010, in American 
Institute of Physics Gonference Series, Vol. 1242, American 
Institute of Physics Gonference Series, ed. G. Bertin, F. de 
Luca, G. Lodato, R. Pozzoli, & M. Rome, 219-230 
Padoan, P., & Nordlund, A. 1997, astro-ph/9706176 
—. 1999, ApJ, 526, 279 

Padoan, P., & Nordlund, A. 2002, ApJ, 576, 870 
—. 2011, ApJ, 730, 40 

Padoan, R, Nordlund, A., & Jones, B. J. T. 1997, MNRAS, 288, 
145 

Pan, L., Padoan, P., Haugbolle, T., & Nordlund, A. 2015, ArXiv 
e-prints 

Pan, L., & Scannapieco, E. 2011, Phys. Rev. E, 83, 045302 
Porter, D., Pouquet, A., Sytine, L, & Woodward, R 1999, 

Physica A Statistical Mechanics and its Applications, 263, 263 


Porter, D., Pouquet, A., & Woodward, P. 2002, Phys. Rev. E, 66, 
026301 

Porter, D. H., Jones, T. W., & Ryu, D. 2015, ArXiv e-prints 
Porter, D. H., Woodward, P. R., & Pouquet, A. 1998, Physics of 
Fluids, 10, 237 

Rathborne, J. M., Johnson, A. M., Jackson, J. M., Shah, R. Y., & 
Simon, R. 2009, ApJS, 182, 131 
Roman-Duval, J., Jackson, J. M., Heyer, M., Johnson, A., 
Rathborne, J., Shah, R., & Simon, R. 2009, ApJ, 699, 1153 
Roman-Duval, J., Jackson, J. M., Heyer, M., Rathborne, J., & 
Simon, R. 2010, ApJ, 723, 492 
Salpeter, E. E. 1955, ApJ, 121, 161 

Sanchez, N., Alfaro, E. J., & Perez, E. 2007, ApJ, 656, 222 
Sanders, D. B., Scoville, N. Z., & Solomon, P. M. 1985, ApJ, 289, 
372 

Schmidt, W., Federrath, G., Hupp, M., Kern, S., & Niemeyer, 

J. G. 2009, A&A, 494, 127 

Schmidt, W., Federrath, G., & Klessen, R. 2008, Physical Review 
Letters, 101, 194505 

Scoville, N. Z., Yun, M. S., Glemens, D. P., Sanders, D. B., & 
Waller, W. H. 1987, ApJS, 63, 821 
Semenov, V. A., Kravtsov, A. V., & Gnedin, N. Y. 2015, ArXiv 
e-prints 

She, Z.-S., & Leveque, E. 1994, Phys. Rev. Lett., 72, 336 
Solomon, P. M., Rivolo, A. R., Barrett, J. W., & Yahil, A. M. 
1987, ApJ, 319, 730 

Stilp, A. M., Dalcanton, J. J., Warren, S. R., Skillman, E., Ott, 

J., & Koribalski, B. 2013, ApJ, 765, 136 
Tamburro, D., Rix, H.-W., Leroy, A. K., Low, M.-M. M., Walter, 
F., Kennicutt, R. G., Brinks, E., & de Blok, W. J. G. 2009, AJ, 
137, 4424 

Teyssier, R. 2002, A&A, 385, 337 

Walch, S., Girichidis, P., Naab, T., Gatto, A., Glover, S. G. O., 
Wiinsch, R., Klessen, R. S., Glark, P. G., Peters, T., Derigs, D., 
& Baczynski, G. 2015, MNRAS, 454, 238 
Walch, S., & Naab, T. 2015, MNRAS, 451, 2757 
Walch, S. K., Girichidis, P., Naab, T., Gatto, A., Glover, 

S. G. O., Wiinsch, R., Klessen, R. S., Glark, P. G., Peters, T., 

& Baczynski, G. 2014, ArXiv e-prints 
Williams, J. R, & McKee, G. F. 1997, ApJ, 476, 166 
Wolflre, M. G., Hollenbach, D., McKee, G. F., Tielens, 

A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152 
Wongwathanarat, A., Muller, E., & Janka, H.-T. 2015, A&A, 577, 
A48 


